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Summary 


The  advection-diffusion  (or  dispersion)  of  gases  from  a  localized,  stationary  or  a  moving  gas 
source  into  an  ambient  environment,  results  in  a  plume  (or  jet)  that  is  representative  of 
processes  with  numerous  applications  of  interest  to  the  Air  Force.  For  example,  the  accidental 
or  deliberate  release  of  gases  from  a  land-based,  air  or  space  vehicle,  results  in  a  plume  that 
can  be  used  for  the  detection  of  the  vehicles  position.  The  deliberate  release  of  biochemical, 
results  in  a  plume  that  can  be  used  in  tracking  and  identifying  the  source.  The  release  of 
gases/odors  from  biological  systems  results  in  a  plume  that  may  be  important  in  search/rescue 
missions.  With  this  work  we  propose  to  further  develop  a  model-based  approach  for  detecting 
the  location  of  a  moving  gaseous  source  and  the  concentration  via  the  use  of  a  Sensing  Aerial 
Vehicle  (SAV).  This  model-based  approach  incorporates  the  dynamics  of  the  advection- 
diffusion  process,  the  concentration  estimation  scheme,  the  sensing  aerial  vehicle  dynamics 
and  guidance,  and  the  onboard  sensor  modeling. 

The  goals  of  research  performed  under  FA9550-12-1-01 14  were  to: 

•  Develop  a  theoretical  and  finite-dimensional  approximation  framework  that  strongly 
couples  theoretical  estimation  and  control  with  advanced  computational  fluid 
dynamics  methods. 

The  objectives  of  research  performed  under  FA9550-12-1-01 14  were  to: 

•  Develop  an  abstract  theoretical  framework  of  the  state  estimator  based  on  a  modified 
Luenberger  observer  with  a  collocated  filter  gain. 

•  Develop  and  implement  in  3D  an  adaptive,  multi-grid,  multi-step  finite  volume 
process-state  estimator.  Upwind  and  flux-limiting  schemes  will  be  incorporated  to 
address  the  convective  and  diffusive  scales  effectively.  The  grid  will  be  constructed 
based  on  length-scales  obtained  through  the  state  estimator  and  adapted  in  regions  of 
interest.  The  time  varying  matrix  representations  of  the  diffusion  operator  due  to  grid 
adaptation  results  in  a  hybrid  dynamical  system. 

•  Incorporate  the  3D  dynamics  of  the  SAV  and  a  sensor  model  with  finite  spatial  and 
temporal  resolution  into  the  SAV  guidance  schemes.  The  Lyapunov  function 
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augmented  by  the  vehicle  dynamics,  links  the  SAV  guidance  to  the  state  estimator 
performance. 


1.  Overview 

The  goal  of  the  work  is  to  provide,  in  real  time,  an  estimate  of  the  gas  concentration 
associated  with  an  emitting  moving  source.  Minimizing  the  damage  from  a  toxic  release  must 
address  not  only  the  source  location  or  its  proximity,  but  also  the  contaminated  material  that 
has  already  been  released  (i.e.  estimate  the  concentration  field).  The  strategy  applied  in  this 
research  effort  is  not  to  reposition  the  sensors  onboard  the  UAVs  to  spatial  areas  of  higher 
concentration  (i.e.  a  local  maximum  concentration),  but  to  send  the  UAVs  with  the  onboard 
concentration  sensors  to  areas  of  higher  state  estimation  error.  An  overview  of  the  proposed 
scheme  is  presented  in  Figure  1. 


Sensing  Unmanned 
Aerial  Vehicles 


Figure  1.  Schematic  of  the  considered  problem  with  an  aerial  source  and  a  group  of  UAVs. 
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2  Physical  Model 


The  mathematical  model  for  the  process  of  gas  release  by  a  moving  point  source  is 
presented  along  with  the  mathematical  model  for  the  sensor  measurement  and  estimator  and 
the  UAV  dynamical  model  and  guidance  law  are  described. 

2.1  Process  Model 

The  plumes  of  interest  in  this  dissertation  are  released  in  the  lower  layers  of  the 
atmosphere,  at  altitudes  below  10  km.  Atmospheric  turbulence  is  responsible  for  the  transfer 
of  the  trace  (plume)  gases  in  the  lowest  layers  of  the  atmosphere  [48],  [49],  [50].  Following 
[24],  [25],  we  present  below  the  “process  model”  which  is  represented  by  the  3D  advection- 
diffusion  equation.  Consider  a  source  of  gas  release,  which  is  moving  along  an  unknown 

trajectory  inside  a  spatial  domain  =  0,  x  hL\  X  [0,Z/^]  G  .  The  spatial  distribution 

of  the  gas  source  is  given  by  the  3D  Dirac  measure  concentrated  at  the  point  of  the  source’s 
spatial  centroid 

b(X,Y,Z)  =  6(X-Xfi))6(Y-  F (t)) 6  (z  -  Zfi)) ,  (2.1) 

where  f  e  is  time.  X,  Y,  and  Z  are  the  spatial  variables,  i^X^(t) ,  F  (t) ,  e  11 .  For 
brevity,  hereinafter  the  time  varying  centroid  of  the  source  is  denoted  by 
0^(f)  =  G  Q.  The  gas  source  is  characterized  by  a  known  release  rate  u(t) 

and  using  Eq.  (2.1),  it  can  be  represented  as 

S{t,e^)  =  b{x,Y,z)u{t).  (2.2) 

Since  the  flows  of  interest  are  turbulent,  the  fluid  velocities  U  ,  V  ,  and  W  in  X ,  Y  , 
and  Z  directions  accordingly  are  random  functions  of  space  and  time.  They  can  be  represented 
as  the  sum  of  deterministic  and  stochastic  component  as  follows 

U  =  U  YU' 

V  =  V  +  V'  (2.3) 

w  ^W  +  W' 

Consider  the  conservation  of  mass  of  gas  within  an  elementary  volume 
17^  =  (ax,  a  Y,  AZj  with  the  center  at  (X,  Y,  zj ,  such  that  0^(i)  G  (7^ ,  shown  in  Figure  3. 
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Assume  the  case  of  chemically  inert  species,  i.e.  the  net  amount  of  material  convected  out  of 
the  volume  element  must  be  balanced  only  by  an  equivalent  amount  of  material  that  is  emitted 

by  the  sources  and  that  enters  by  molecular  diffusion.  The  conservation  of  mass  or 

concentration  c  through  this  elementary  volume  is  written  in  integral  form  as 


!///<='>«  =-^(F-)d^  +  ///5da.  (2.4) 

^  Q  A  Q 

C  C 

where  A  is  the  volume  surface  area,  n  is  the  normal  vector  to  the  surface,  and  F  is  the  total 
flux  through  the  surface.  Applying  the  Ostrogradsky-Gauss  theorem  yields 

=  (2.5) 

A  Q 

C 


The  gradient  operator  in  Cartesian  reference  frame  is  V  ■  F  =  — —  H - —  H - — .  The 

dX  dY  dZ 

flux  through  each  face  F^,  ,  and  F^  of  the  elementary  volume  is  represented  as  a  sum  of 

advective  and  diffusive  fluxes.  According  to  the  Pick's  law  [49],  the  diffusive  flux  is 
proportional  to  the  concentration  gradient.  Therefore,  the  total  flux  through  each  face  is  written 
as  follows 
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—  Oc  —  dc 

F=cU  ax-D—  -  cU  ax-D  — 

^  dX  dX  ^  ^ 

.  2/1  ^2 


F  =  cF  AK  -D^ 

.  2 


-  cy  ^.-D- 


2 


—  dc  —  dc 

F  =  cW  az-D—  -  cW  az-D  — 

^  dz  dz 

V  2  J  V  2  ^ 

where  D  is  molecular  diffusivity.  In  order  to  express  the  fluxes  through  the  faces  X  H 

AY  AZ  AX  AY 

Y  H - ,  and  Z  H - in  terms  of  the  fluxes  through  X - ,  Y - ,  and  Z 


we  apply  the  Taylor  series  expansion 


d{cU) 


dX 

2 


cy  Av=cy 


d{cV)\ 


BY 

cir|  A,  =  <;ir|  +  AZ 

dc  dc  d  (  dc) 

D—  =D—  F—  D—  AX 
dX  ^  ^  dX  ^  ^  dx[  dXj  ^ 

2  2  ^2 

+^L^]  aF 

9F,,^^  ar^  ay),,  „ 

2  2  2 

dc  dc  d  (  dc  1 

dYX  =D—  +—  D—  AZ 

dZ  dZ  dz(  dzj  ^ 

2  2  2 

Eq.  (2.7)  implies  that  the  gradient  of  the  total  flux  is 

V-F  =  I  ^  ^  ^  {d^^\  (2  8) 

dX  dY  dZ  dX  dX  dY  dY  dZ  dZ 

V  /  \  J  \  J 

Substituting  this  result  into  Eq.  (2.4)  supplemented  with  Eq.  (2.5)  and  Eq.  (2.3),  integrating 
and  dividing  by  the  elementary  volume  yields 
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(2.9) 


dc 

dt 


+ 


a  (c  ([/  +  [/'))  a(c((/  +  F'))  dl^c(w  +  w') 


d 


dX 


D 


dX, 

dc 


dX 


+ 


d 


dY 


D 


dY 
dc 


dZ 


dY 


+ 


A 

dZ 


D 


dc 


+  <S(t,0j 


Since  U' ,  V' ,  and  W'  are  random  variables,  concentrations  obtained  from  solution  of  Eq. 
(2.9)  are  also  random  variables.  Thus,  the  determination  of  a  concentration  c  in  the  sense  of  a 
specified  function  of  space  and  time  is  not  possible.  Instead,  the  mean  concentration  (c)  should 


be  specified,  such  that  c  =  (c)  -\-  c  ,  with  (c)  =  0 .  Averaging  Eq.  (2.9)  over  an  infinite 
ensemble  of  realizations  of  the  turbulence  yields 


A) 

dt 


d 

dX 


(U{c)  +  (UV)}  +  ;^iv(c}  +  (W)}  +  ^{W(c}  +  (WV)} 

d 


d 


d 


dx 


D 


d{c) 


dX 


+ 


d 


dY 


D 


d{c) 


dY 


+ 


dZ 


D 


d{c) 


dZ 


(2.10) 


In  order  to  relate  the  turbulent  fluxes  {U'c') ,  {V'c') ,  and  {W'c')  to  (c)  we  apply  the 
mixing-length  model  [48].  The  mixing  length  is  a  measure  of  the  maximum  distance  in  the 
fluid  over  which  the  velocity  fluctuations  are  correlated  (i.e.  a  measure  of  the  eddy  size). 
According  to  the  mixing-length  model. 


(C/'c') 


-K 


XX 


d{c) 

dX 


(E'c')  = 
{W'c')  = 


9{c) 

dY 

d{c) 

dZ 


(2.11) 


where  =  1?2, 3  is  the  eddy  diffusivity  tensor.  In  the  case  considered  here  the 


coordinate  axes  coincide  with  the  principal  axes  of  the  tensor,  therefore  it  has  only  three 
nonzero  diagonal  elements  if .  Assume,  in  addition  that  molecular  diffusion  is 


negligible  compared  with  the  turbulent  diffusion 
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(2.12) 


d 


dX 

d 


D 


d{c) 


dY 

A 

dZ 


D 


dX 

d{c) 


«— ([/V) 


D 


dY 

d{c) 


dZ 


<s; 


dX 

A 

dY 

A 

dX 


(FV) 

(W'c) 


Substituting  Eq.  (2.11)  into  Eq.  (2.10)  with  assumption  (2.12)  yields  the  advection- 
diffusion  equation  written  in  the  so  called  conservative  form. 


a{c)  9(u{c))  a(y(c))  d[w(c)) 


dt 

d 


dX 


•  + 


K 


dX 

d{c) 


'  + 


dY 


+  ■ 


dZ 


XX 


dX 


+ 


d 


dY 


d{c) 


+ 


A 

dZ 


K 


d{c) 


zz 


dZ 


+ 


5(i.©) 


(2.13) 


rr  QY 

which  describes  the  dispersion  of  the  plume  created  by  the  point  gas  source  in  the  ambient 
atmosphere. 

a(c)(t)| 


Equation  (2.13)  is  supplemented  with  Neumann  boundary  conditions 


dn 


=  0  at 


on 


the  boundaries  of  the  domain  Q  and  an  initial  condition  {c){t  —  0) 


=  0. 


In  case  of  constant  wind  speeds  U ,  V ,  W  and  eddy  diffusivities  the 


advection-diffusion  equation  (2.13)  is  written  in  a  strong  conservative  form  as 


dt 


dX  dY  dZ 
d\c)  d\c)  d\c) 

dX^  dY^  dZ^ 


+  5(t,  ©J 


(2.14) 


In  this  work,  the  atmospheric  wind  speeds  and  eddy  diffusivities  are  assumed  to  be  known 
functions  of  the  spatial  variables.  As  was  mentioned  above,  the  process  of  gas  release  is 
considered  in  the  lowest  layer  of  the  atmosphere,  called  the  planetary  boundary  layer  (PEL). 
The  thickness  of  the  PEL  is  a  variable  in  both  time  and  space  and  varies  from  20-500m  during 
nighttime  hours  to  0.2-5  km  in  the  late  afternoon.  Since  the  time  scale  and  space  scale  of  the 
process  under  consideration  is  small  (up  to  10-15  min  and  up  to  4  km),  the  thickness  of  the 
PEL  in  this  case  is  considered  to  be  constant  and  equal  to  0.5-1  km.  The  mean  wind  speed  in 
the  PEL  generally  increases  with  height.  An  empirical  representation  of  the  mean  wind 
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distribution,  which  is  frequently  used  in  air  pollution  dispersion  applications  is  the  power-law 
profile 


U 

r 


{  \m 

z] 


m  <1 


(2.15) 


where  is  the  wind  speed  at  the  reference  height  and  m  is  an  exponent  less  or  equal  to 

unity.  The  parameter  m  depends  on  the  surface  roughness  and  stability.  Observations  have 
shown  that  surface  roughness  and  stability  increase  when  m  increases.  Under  near-neutral 
conditions  m  ranges  from  0.15  for  smooth  water,  snow  and  ice  surfaces  to  0.4  for  urban  areas. 
Eddy  diffusivities  and  are  usually  either  assumed  to  be  constants  or  specified 

as  functions  of  height  through  power-law  relations  of  the  type 


K 


r 


V  r ; 


n  <1 


where  parameter  n  depends  on  the  surface  roughness  and  stability. 


(2.16) 


2.2  Sensor  Model 

Various  sensing  techniques  exist  to  measure  concentration  of  trace  gases.  These  principles 
are  based  on  the  physical  and  chemical  properties  of  the  gases.  Basic  sensor  characteristics  are 
sensitivity,  range,  precision,  accuracy,  resolution,  response  time,  offset,  and  hysteresis.  Among 
sensing  technologies,  the  most  commonly  used  are  electrochemical,  solidelectrolite,  catalytic, 
spectroscopic  and  photoionization  sensors  [51],  [52]. 

In  this  work,  the  detection  and  quantification  of  pollutants  and  natural  trace 
environmental  chemicals  are  of  interest.  This  requires  sensors  with  high  level  of  sensitivity. 
Furthermore,  the  required  sensor  should  have  very  small  response  time  (within  a  couple  of 
seconds)  to  provide  measurements  for  the  estimation  process.  This  is  due  to  the  fact  that  the 
guidance  scheme  is  coupled  to  the  estimator  performance  and  the  resulting  control  inputs  for 
the  UAV  should  be  processed  in  real  time.  Therefore,  spectroscopic  sensing  technology  seems 
to  be  the  most  acceptable  for  this  work. 

The  principle  of  operation  of  a  spectroscopic  sensor  is  based  on  a  unique  absorption 
spectrum  exhibited  by  each  molecule.  The  spectrum  is  described  in  terms  of  absorption 
strength  versus  frequency,  temperature  and  pressure  in  spectroscopic  databases.  The  molecular 
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structure  determines  at  which  frequencies  the  incident  light  will  be  absorbed  by  the  molecule. 
Among  optical  sensing  techniques  the  infrared  (IR)-source  gas  sensors  are  widely  used. 

An  IR-source  gas  sensor  contains  three  major  parts,  shown  in  Figure  4:  the  IR  source, 
a  volume  of  the  gas  sample  and  the  IR  detector.  When  the  IR  source  emits  broadband  radiation 
including  the  wavelength  absorbed  by  the  target  gas,  the  sample  gas  in  the  gas  cell  will  absorb 
the  radiation  in  its  particular  way.  The  optical  filter  is  used  to  screen  out  all  radiation  except 
for  the  wavelength  that  is  absorbed  by  the  target  gas.  Therefore,  the  presence  of  interested  gas 
could  be  detected  and  measured  by  an  IR  detector.  This  system  is  also  known  as  Non- 
Dispersive  Infrared  (NDIR)  gas  sensor  [52].  The  smallest  NDIR  sensor  modules  weight 
approximately  5-10  grams  and  have  size  less  than  40  mmx40  mmxl5  mm. 

The  estimation  process  is  based  on  in  situ  sensing,  i.e.  a  sensing  technique  where  a 
device  is  in  direct  contact  with  the  environmental  phenomena  (here,  concentration).  In  situ 
sensing  in  fluid  environments  is  classified  into  Eulerian  and  Lagrangian  based  on  the  different 
reference  frames  [53].  In  the  present  work,  the  Lagrangian  technique  is  used,  according  to 
which  sensors  move  freely  in  the  fluid  itself,  and  gather  measurements  as  they  moves  through 
the  environment  [54],  [55]. 

Mathematically,  the  gas  sensor  location  0^  G  f)  is  described 

similar  to  the  gas  source  location  (see  Eq.  (2.1))  with  the  3D  Dirac  delta  function.  Therefore, 
the  sensor  readings  is  provided  via 

Ly  Ly 

y[t,e\t))  =  f  f  f(c}(t,X,Y,z}6(X-X^(t))s(Y-V(t))s(Z-Z(t))dZdYdX  (2.17) 

0  0  0 

Sample  Gas  Input  Sample  Gas  Output 


Detector 


IR  Source  Gas  Chamber  Optical  Filter 

Figure  4.  IR-source  gas  sensor  based  on  the  basic  absorption  spectrometry  [52].  Copyright  2012  by 
the  authors;  licensee  MDPI,  Basel,  Switzerland 
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2.3  Estimator  Model 


The  estimator  model  is  based  on  the  advection-diffusion  equation  (2.13).  It  takes  the  form 
of  a  Luenberger  observer  supplemented  with  an  output  injection  term  72.(t,0  )  that  is 

dependent  on  the  sensor  location  : 


dt  ^ 

dX 

-  +  — 

d 

0(6)] 

dX 

dX 

dY 

+ 


d 


+ 

( 


dY 


dZ 

d{c) 


+ 


A 

dZ 


K 


d{c) 


zz 


dZ 


(2.18) 


+  7l{t,Q 


YY  QY 

where  (c)  is  the  estimated  mean  concentration.  The  boundary  conditions  for  Eq.  (2.18)  are 

d{c){t)\ 


dn 


0  and  the  initial  conditions  are  {c){t  =  0)  =  {c)  {t  =  0)  .  The  output  injection 


an 


term  1Z{t,  0  J  is  specified  by  the  difference  between  a  “true”  concentration  and  a  state  estimate 

at  the  current  UAV  location,  multiplied  by  the  filter  gain.  The  latter  is  taken  to  be  a  weighted 
multiple  of  the  dual  of  the  observation  operator  associated  with  the  sensor's  spatial  distribution 
given  in  Eq.  (2.17) 

7^  (t,  e^(t))  =  T6(X-  X  (t))  6(Y-Y^ (t))  S(Z-Z^ (t))  [y  (t,  0  )  -  (c)  (t,  0^ )]  (2. 19) 
where  T  >  0  is  a  user-defined  estimation  gain. 


2.4  UAV  Dynamics  and  Guidance 


2.4.1  Review  ofUAVs 

Unmanned  Aircraft  Vehicles  (UAVs)  have  evolved  rapidly  over  the  past  decade  driven 
primarily  by  military  uses,  and  have  begun  finding  application  among  civilian  users  for  earth 
sensing  reconnaissance  and  scientific  data  collection  purposes  [43].  Among  advantages  of  the 
UAVs  as  compared  to  manned  aircrafts  is  long  flight  duration,  improved  mission  safety,  flight 
repeatability  due  to  improving  autopilots,  and  reduced  operational  costs. 


Among  the  low-cost  and  low-risk  UAVs  for  remote  sensing  expeditions,  one  should  mention 
the  Aerosonde®  UAV.  This  UAV  can  be  programmed  to  make  very  detailed  flight  patterns  that 
can  be  flown  automatically  and  in  very  extreme  weather  conditions.  Specifications  for  the 
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Aerosonde®  UAV  are  provided  below,  as  well  as  a  discussion  of  its  applications  and 
advantages  [45]. 

2.4.2  Coordinate  frames 

To  derive  the  dynamic  model  for  the  UAV,  the  inertial  and  the  body  coordinate  systems 
(or  frames)  are  required  [56],  [57],  [58].  The  inertia  frame  is  used  in  applying  Newton’s  law 
and  is  the  frame  in  which  GPS  provides  position  and  speed.  The  body-fixed  frame  is  often  used 
for  describing  the  aerodynamics  forces  as  well  as  for  on-board  sensors,  such  as  accelerometers 
and  rate  gyros.  Therefore  it  is  important  to  identify  these  coordinate  frames  and  to  describe  the 
transformation  between  them. 

For  a  short-range  UAV  (with  the  flight  duration  of  up  to  10  hours  and  maximum  distance 
of  several  kilometers)  the  earth  is  assumed  to  be  flat  and  non-rotating.  The  inertial  coordinate 
frame  is  an  earth-fixed  coordinate  system  denoted  byjT*  and  is  depicted  in  Figure  5.  The  unit 
vectors  X,  Y,  and  Z  are  directed  north,  east  and  towards  the  earth  center  respectively.  The 

body  coordinate  frame,  which  is  denoted  hy  ,  has  the  origin  at  the  center  of  mass  of  the 
vehicle.  The  axes  are  directed  as  follows;  x  points  out  the  nose  of  the  airframe,  y  points  out 

the  right  wing,  and  z  points  out  the  belly.  The  body  frame  is  obtained  from  the  inertial 
frame  T'  by  the  four  consecutive  transformations,  depicted  in  Figure  5. 

•  Translation  of  the  origin  of  the  T"  reference  frame  to  the  center  of  mass  of  the  vehicle 
to  result  in  the  vehicle  frame  T" 

•  Positive  right-handed  rotation  about  the  z”  axis  by  the  heading  (yaw)  angle  ip  to  result 
in  the  vehicle- 1  reference  frame 

•  Positive  right-handed  rotation  about  the  y*'^  axis  by  the  pitch  angle  9  to  result  in  the 
vehicle-2  reference  frame 

•  Positive  right-handed  rotation  about  the  axis  by  the  roll  (bank)  angle  (p . 
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Figure  5.  Transformation  from  the  inertial  frame  to  the  body  frame. 


The  transformation  of  an  arbitrary  vector  p  from  the  vehicle  frame  T"  to  the  body  frame 


is  given  by 


R*  y,,  0,  V-)  p-  =  K  M  K  (»)  K  M  p 
0 


1  0 
0  cos  (p  sin  (p 

0  —  sin  (p  cos  (p 


cos  ^  0  —  sin  9 
0  1  0 

sin  0  0  cos  0 


cos  Ip  sin  pj  0 
—  sin  Ip  cos  Ip  0 


0 


C^C,  c^s,  -s^ 

9  Ip  9  Ip  9 


+  C^S^S^-Sf^ 


0 


(2.20) 


where  'RP  {cp^  9,  r/)  j  is  the  transformation  matrix  from  the  vehicle  frame  to  the  body  frame,  71^ 


is  the  transformation  matrix  from  the  reference  frame  to  the  body  reference  frame, 
is  the  transformation  matrix  from  the  reference  frame  to  the  reference  frame  ,  7^”^  is 

’  V 


the  transformation  matrix  from  the  vehicle  reference  frame  to  the  reference  frame  , 
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X 


Figure  6.  Transformation  from  the  body  frame  to  the  wind  frame. 


C  =  cos*  and  S  =  sin* .  The  matrix  in  Eq.  (2.20)  represents  Euler  transformation  matrix 

>):  >):  'D  ^  ^ 

with  Euler  angles  0,  0,  and  V'.  The  rotation  sequence  (f)-9-'ip  is  commonly  used  for 
representing  aircraft  orientation  in  three  dimensions. 

The  wind  direction  relative  to  the  body  reference  frame  is  specified  using  the  angle  of 
attack  a  and  the  sideslip  angle  (3 .  The  angle  of  attack  a  is  defined  as  a  left-handed  rotation 

about  the  body  y  axis  to  result  in  the  stability  reference  frame  T'  as  shown  in  Figure  6  (left). 
The  sideslip  angle  (3  then  results  from  the  right-handed  rotation  of  the  stability  frame 


about  z®  to  yield  the  wind  frame  JT™  as  shown  in  Figure  6  (right).  The  unit  vector  x  of  the 
wind  frame  JT®  is  aligned  with  the  airspeed  vector  which  is  the  velocity  of  the  airframe 
relative  to  the  surrounding  air.  Therefore,  the  transformation  of  an  arbitrary  vector  p  from  the 


body  frame  to  the  wind  frame  JT™  is  given  by 


)p‘ 

II 

a)n:{l3) 

COS  a 

0 

sin  a 

COS  [3 

sin  [3 

o' 

0 

1 

0 

—  sin  (3 

COS  (3 

0 

—  sin  a 

0 

cos  a 

/ 

0 

0 

1 

/ 

with 


cos  f3  cos  a  sin  (3 
—  sin  (3  cos  a  cos  (3 
—  sin  a  0 


cos  f3  sin  a 
—  sin  (3  sin  a 
cos  a 


(2.21) 


(2.22) 
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2.4.3  Wind  Triangle 

For  a  small  UAV,  whose  speed  ranges  from  20  to  35  m/s,  the  wind  speed  ranges  from  20 
to  50  percent  of  the  airspeed,  therefore,  it  is  important  to  account  for  the  effect  of  wind  on  the 
UAV  dynamics.  The  inertial  forces  acting  on  the  UAV  are  dependent  on  the  velocity  relative 
to  a  fixed  inertial  reference  frame  T'  referred  to  as  ground  speed  V^ ,  whereas  the 

aerodynamic  forces  depend  on  the  airspeed  V^ .  Denoting  the  wind  velocity  relative  to  the 
inertial  frame  by  V^ ,  these  velocities  are  related  by  the  expression 


V  =  V  -  V  (2.23) 

The  relationship  (2.23)  is  called  the  wind  triangle.  Given  the  wind  speed  components 
U,  V,  and  W  with  respect  to  the  inertial  reference  frame,  i.e.  north,  east,  and  down,  and  using 


Eq.  (2.20),  the  wind  speed  can  be  expressed  in  the  body  frame  JT''  as 


V^  = 


/  \ 

u 

/  \ 

U 

w 

V 

V 

w 

w 

V  J 

W 

[  w) 

(2.24) 


Defining  u,  v,  and  w  as  the  body  frame  components  of  the  ground  speed  V^ ,  and  u^,  v^, 


and  as  the  body  frame  components  of  the  airspeed  V^ ,  Eq.  (2.23)  is  written  as 


/  \  /  \ 
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u  —  U 
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w 

y’’  = 
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— 

V  —  V 
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w 
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w  —  w 

(2.25) 


On  the  other  hand,  the  airspeed  vector  in  the  reference  frame  related  to  the  wind  JT®  has 
only  one  nonzero  component 


(2.26) 


Therefore,  combining  Eq.  (2.22)  and  Eq.  (2.25)  yields 
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Inverting  Eq.  (2.27)  gives 


V  =  Ju^  +  r;'  +  w;' 

a  yj  r  r  r 


a 


/5 


tan 


-1 


^  \ 
w 


u 


sin 


-1 


2  I  2  I  2 

M  +  r;  +  w 


(2.27) 


(2.28) 


The  direction  of  the  ground  speed  vector  relative  to  the  inertial  frame  T'  is  specified  using 
the  course  angle  x  and  the  flight  path  angle  7  (see  Figure  7).  The  course  angle  x  is  defined 


as  the  angle  between  the  horizontal  component  of  the  ground  speed  vector  and  the  vector 
X  (true  north);  the  flight  path  angle  7  is  defined  as  the  angle  between  the  ground  speed 
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and  its  horizontal  component.  The  ground  speed  vector  in  the  inertial  frame  can  be  expressed 
in  terms  of  %  and  7  as 


V* 

9 


COS  X  cos  7 
sin  X  cos  7 
—  sin  7 


(2.29) 


The  geometrical  representation  of  the  concept  of  wind  triangle  is  given  in  Figure  8.  Here 
X^  is  referred  to  as  crab  angle  and  7^  is  the  air-mass-referenced  flight-path  angle: 


x.-x-i’ 

7  =e-a 

'  a 


(2.30) 


Assuming  sideslip  angle  (3  to  be  negligible  so  that  the  airspeed  vector  is  aligned  with 


the  X  axis  (Figure  8),  the  airspeed  can  be  expressed  in  the  inertial  frame  as 


Y'  =V 


cos  'Ip  cos  7^ 
sin  'ip  cos  7^ 
—  sin  7 

I  a 


(2.31) 


Therefore,  the  wind  triangle  (2.23)  can  be  rewritten  as 


/  \ 

COS  X  cos  7 
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cos  Ip  cos  7^ 
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sin  X  cos  7 
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=  V 

a 

sin  Ip  cos  7^ 

—  sin  7 
\  J 

w 

\  J 

—  sin  7 

(2.32) 


Taking  squared  norm  of  both  sides  of  Eq.  (2.32)  yields 
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(2.33) 
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sin  X  cos  7 
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The  third  row  of  Eq.  (2.32)  is  solved  for  the  air-mass-referenced  flight  path  angle  7 


7  =  sin 

'  a 


V  Sin7  +  1F^ 

n  ' 


V 


(2.34) 


Multiplying  both  sides  of  Eq.  (2.32)  by  sin  %,  cos  %,  0 j  and  solving  for  the  heading 

angle  1/^  yields 


V’  =  X  -  sin  ^ 
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(  \ 
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-  sin  X  j 

V  cos  7 

a  '  a 
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\  J 

cosx 

V  / 

/ 

(2.35) 


2.4.4  Dynamic  Model 

The  UAV  equations  of  motion  are  derived  based  on  commonly  encountered  in  the 
literature  point-mass  model  of  a  fixed-wing  aircraft  [56],  [60],  [61],  [62].  This  model  is  drawn 
from  free-body  diagrams  and  includes  lift,  drag,  and  thrust  forces. 

Figure  9  provides  free-body  diagrams  for  a  UAV  in  climbing  coordinated  turn:  to  the  left 
the  UAV  is  rolling  at  an  angle  of  </>  and  forces  are  shown  in  x  —  z  plane,  and  to  the  right  the 
diagram  shows  the  view  in  the  direction  of  — x  axis,  forces  are  shown  in  y  -  z  plane.  Thrust 
force  is  denoted  by  T ,  drag  and  lift  forces  are  denoted  by  D  and  L  respectively;  M  is  the 
mass  of  the  UAV.  Applying  Newton's  second  law  along  the  x  axis  (Figure  9,  left)  and  in  the 
vertical  direction  (Figure  9,  right)  yields 


Figure  9.  Free-body  diagram  of  the  forces  acting  on  the  UAV  in  the  climbing  coordinated  turn. 
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—  T  —  D  —  Mg  sin  7 
Mg  cos  7  —  L  cos  0  =  0 


(2.36) 

(2.37) 


During  a  coordinated  turn,  the  bank  angle  0  is  set  so  that  there  is  no  net  side  force  on  the 
UAV.  In  order  to  derive  relationship  for  course  angle  x  terms  of  lift  L,  eonsider  the 
coordinated  turn  maneuver  from  the  top  view  of  horizontal  plane,  shown  in  Figure  10.  The 
horizontal  eomponent  of  the  lift  foree  L  is  acting  in  the  radial  direction  opposite  to  the 
centrifugal  force,  therefore  the  Newton's  second  law  results  in 


M  cos  7  j  X  ~  sin  0  cos  [x  ~  V')  =  0 


(2.38) 


Solving  Eq.  (2.38)  for  x  provides 


L  sin  0  cos  (x  ~  0) 
MV  cos  7 

9  ' 


(2.39) 


To  derive  the  dynamics  for  the  flight  path  angle  7 ,  consider  a  pull-up  maneuver  in  which 
the  aircraft  climbs  along  an  arc  as  shown  in  Figure  10,  so  that  the  airframe  is  rolled  at  an  angle 
0 .  Newton's  second  law  in  the  radial  direetion  gives 


L  cos  0  —  Mg  cos  7  —  MV  7  =  0 


(2.40) 


Solving  Eq.  (2.40)  for  7  gives 


L  ,9 

7  = - cos© - cos  7 

MV  V 

9  9 


(2.41) 


Combining  Eq.  (2.36),  Eq.  (2.39),  and  Eq.  (2.41)  with  the  kinematie  equations  (2.29) 
yields  the  equations  of  motion  that  describe  the  point-mass  model  of  the  fixed-wing  aireraft: 

X  =  E  cos  X  cos  7 
Y  =  V^  sin  X  cos  7 
Z  =  —V  sin  7 

V =^-^-  ■ 

^  M  M  (2.42) 

L  ±9 

7  = - cos  © - cos  7 

MV  V 

9  9 

L  sin  0  cos  i^x  ~  0) 

^  MV  cos  7 

a  ' 


with  lift  and  drag  being  determined  as 
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(2.43) 


L  =  -pV^SC^ 

2  “  ^ 

D  =  -pV^SC^ 

2  “  ^ 

where  (7^  and  (7^  are  lift  and  drag  coefficients  related  as  follows 


C 


D 


+ 


ireAR 


(2.44) 


A  0 

where  C  is  the  parasitic  drag  due  to  the  shear  stress  of  air  moving  over  the  wing;  AR  =  — 
”  S 

is  the  wing  aspect  ratio;  b  is  the  wingspan,  and  e  is  the  Oswald  efficiency  factor  [56],  [58]. 

The  lift  coefficient  for  the  point-mass  model  (2.42)  is  given  by  the  linear  function  of  the  angle 

of  attack: 

(a)  =  +C^a  (2.45) 


where  C  is  the  lift  coefficient  at  zero  angle  of  attack  and  C  = - -  is  the  stability 

»  “  da 

derivative.  The  magnitudes  for  C  and  C  are  found  through  wind  tunnel  tests  or  a  detailed 

u  a 

computational  study  and  usually  can  be  found  in  the  aircraft  specifications. 


Figure  10  Free-body  diagram  of  the  forces  acting  on  the  UAV  in  the  climbing  coordinated  turn:  top 
view  of  the  horizontal  plane  and  side  view  of  the  vertical  plane. 
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The  engine  thrust  T ,  the  bank  angle  cf) ,  and  the  lift  coefficient  C ^  are  the  control  inputs 


to  the  dynamic  model  (2.42).  The  control  inputs  are  assumed  to  stay  within  certain  limits 


0  <  r  <  r 

—  —  max 


< 


0< 


(2.46) 


where  T  ,  6  ,  and  C,  are  pre-specified  constants. 

The  model  (2.42)  is  popular  for  simulation  of  the  UAV  motion  [60],  [61],  [62]  due  to  the 
fact  that  it  models  the  aircraft  response  to  inputs  that  a  pilot  commonly  controls:  engine  thrust 
T ,  lift  coefficient  C ^ ,  and  bank  angle  0 .  At  the  same  time,  it  allows  to  avoid  cumbersome 

relations  in  the  12-state  model  for  the  UAV  kinematics  and  dynamics,  which  is  usually  used 
to  implement  autopilot  for  specific  aircraft  maneuvers. 


2.4.5  UAV  Guidance 


In  order  to  derive  the  UAV  guidance  scheme,  it  is  convenient  to  rewrite  the  process  model 
equation  (2.13),  and  the  estimator  equation  (2.18)  in  the  abstract  form  as  follows 

x^=^  +  ,S(f,0j  (2.47) 

=  .4x -|- (^*rc  —  X  j  (2.48) 

where  x^fj  =  is  the  concentration  state,  x^^j  is  the  estimated  concentration  state. 


A  is  the  advection  diffusion  operator,  which  is  specified  based  on  Eq.  (2.13) 


0 

ivu) 

8  A) 

d[p>w) 

dX 

dY 

dZ 

d 

dx 

dX 

) 

dY 

dY 

dZ 

dZ 

.  / 

(2.49) 


for  (p  &  and  C  is  the  output  operator  associated  with  the  sensor's  spatial  location 

Lx  Ly  Ly 

Cip^fff  cpS(x  -  X^(t))s(Y  -Y(t))s(z  -  Z^(t))dZdYdX  (2.50) 

0  0  0 

for  all  the  test  functions  (p  &  .  The  UAV  guidance  scheme  is  based  on  the  estimator 

performance  through  the  state-estimation  error 


20 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


eit,X,Y,z'^  =  {c){t,X,Y,Z^  —  {c)it,X,Y,z'^ .  Combining  Eq.  (2.47)  and  Eq.  (2.48)  yields 
the  governing  equation  for  the  evolution  of  the  state  estimation  error 

=  (^  -  c*rc) e  +  S(t,e^)  =  A^e  +  S{t,e^)  (2.5 1) 


Equation  (2.51)  is  supplemented  with  Neumann  boundary  conditions 


initial  conditions  are  e{t  =  0)  =  e  {t  =  0) 


de{t) 


dn 


—  0 ,  and  the 


an 


The  control  inputs  should  be  chosen  so  that  the  UAV  is  driven  towards  the  spatial  areas 
of  higher  estimation  error.  For  this  purpose,  the  Lyapunov  redesign  method  is  used  with  the 
following  choice  of  the  Lyapunov  functional  [63] 

e  =  -{e,\e)  (2.52) 

where  (•,  •)  denotes  the  (o)  inner  product.  The  system  is  stable  provided  that  the  derivative 

of  the  Lyapunov  functional  in  Eq.  (2.52)  is  negative  semidefinite,  i.e.  5  <  0 .  Taking  time 
derivative  of  Eq.  (2.52)  yields 


£^-2 


Axh)  =  -MZ  +  +  ^^4  +  M)  (2-53) 


dt 


where 


=  e,  (f)  = 


^d{c){t,e 


d 


d 


*  _ 


^X,Y,Z.  (2.54) 


To  ensure  the  derivative  of  the  Lyapunov  functional  (2.53)  is  negative  semidefinite,  the 
UAV  must  have  the  following  Cartesian  velocities 

=  -Kee^ 

Y^  =  -k^ee^  (2.55) 

Z'*  =  -Kee^ 

where  >  0,  ky>Q ,  and  k^>  Q  are  user-defined  constant  guidance  gains. 

The  desired  Cartesian  components  of  the  velocity  vector  in  Eq.  (2.55)  yield  the  desired 
magnitude  of  the  ground  speed  U  ,  the  course  angle  %  and  the  flight  path  angle  7  are 

calculated  using  the  second  and  the  third  equations  of  the  dynamic  model  (2.42)  accordingly 
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= 


+  Y‘‘^  +  Z 


d2 


7 


arcsin 


X*"  =  arctan  2 


Zd 

yd 

9 

Yd 


(2.56) 


Xd 


Given  the  values  of  desired  ground  speed  and  current  ground  speed,  the  desired 
acceleration  is  defined  as  a  numerical  approximation  of  the  time  derivative;  the  flight  path  rate 
and  the  course  rate  are  represented  by  the  first-order  models  [56] 


• .  V^-V 

yd  _  g  g 

^  At 

Y  =  [Y  -  Yj 

Y 


(2.57) 


K  [Y  -  x) 

where,  >  0  and  >  0  are  user-defined  constants.  Substituting  Eq.  (2.57)  into  the  dynamic 
equations  (2.42)  supplemented  with  Eq.  (2.43),  yields  the  control  inputs  to  the  UAV 

0  =  arctan  2 


f 

/  \ 

..  1  ffcosY 

'  U'* 

\ 

Y  cosY, 

cos(Y  -  Y) 

,  g  j 

MV‘‘ 


c,= 


'  v‘‘ 


(2.58) 


X  pV^S  cos  (P 
2 


T  =  MV^  +  -  pV^S  ( +  KCl )  +  Mg  sin  7“ 
2 


The  UAV  is  subject  to  realistic  constraints  and  therefore,  the  control  inputs  (2.58)  should 
stay  within  certain  limits  (2.46).  The  linear  aerodynamic  model  is  chosen  for  the  lift  and  drag, 

therefore,  67““  is  set  to  be  equal  to  the  lift  coefficient  at  stall  i.e., 

67““  =  (2.59) 
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The  maximum  thrust,  is  provided  by  the  performance  characteristics  of  the  particular 
UAV  used.  In  order  to  obtain  the  maximum  bank  angle,  d)  ,  the  load  factor  defined  as  the 

^  “max 

ratio  of  the  lift  to  the  weight  of  the  aircraft  is  used 


n 


L 

Mg 


(2.60) 


For  a  small  UAV  in  turning  level  flight,  the  maximum  load  factor  is  usually  equal  to 
1.5  [59].  Therefore,  the  condition  for  results  from  the  force  balance  equation  (fifth 
equation  in  system  (7))  with  7  =  7  =  0 

1  L 


COS0 


(2.61) 


2.4.6  Extension  to  Multiple  Sensors 


The  desired  Cartesian  velocities  in  Eq.  (2.55)  are  expressed  in  terms  of  the  estimation  error 
£  and  the  error  gradients  e*  at  the  current  sensor  location  0^  .  Therefore,  in  order  to 


implement  the  proposed  guidance  scheme,  we  should  have  at  least  seven  point  sensor 
measurements.  These  sensors  can  be  either  attached  to  a  single  UAV  or  to  different  UAVs  that 
maintain  a  rigid  flying  formation.  The  gradient  is  then  calculated  using  a  central  difference 
approximation  of  the  spatial  derivative  based  on  the  six  additional  measurements  (two  in  each 
Cartesian  direction  X ,  Y  ,  and  Z ). 

It  seems  reasonable  to  include  additional  measurements  into  the  estimator  (2.18).  For  this 
purpose,  the  observation  operator  C  in  Eq.  (2.48)  is  given  by  the  7-dimensional  vector 


Cif  — 


0  0  0 

Lx  Ly  Ly 

[000 


Ly  Ly  Ly 

III 


X  (f))  5  (f  -  F  (O)  ^  [Z  -  Z  (f))  dZd  FdX 


X ))  (5  (y  -  F  (t ))  6  (f  -  Z  ^  (t ))  dZd  FdX 


(2.62) 


for  all  ip  E  (d).  The  filter  gain  matrix 
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Figure  11.  Example  of  the  flying  formation  of  UAVs  for  the  concentration  gradient  measurements. 


r 


rn  0  0 

0  0 

0  0  r„ 


(2.63) 


and  which  is  taken  to  be  diagonal  for  simplicity.  A  full  matrix  would  designate  an  all-to-all 
communication  between  the  7  UAVs,  require  7  interacting  estimators.  The  diagonal  structure 
of  the  filter  gain  matrix  F  ensures  that  there  is  a  single  centralized  filter  with  the  6  follower 
UAVs  in  strict  formation  about  the  leader  UAV  (marked  as  s4).  Therefore,  the  derivative  of 
the  Lyapunov  functional  (2.53)  yields 


f=-2|.4„ef  +  r.( 


EE^X.  +  EEyY  +  EE^Z^ 


+  ...  +  (^EEj^X  +  EEyY  +  EE^Z^ 


si 


(2.64) 


where  the  error,  the  error  gradients,  and  the  Cartesian  velocities  are  calculated  at  each  sensor 
location  denoted  by  si, ...  57 . 

Consider  a  flying  formation,  where  all  the  vehicles  maintain  a  constant  distance  from  the 
leader  UAV  (denoted  by  s4),  shown  in  Figure  11.  In  this  case 

X  =  V 

SI  sj 

Y,=y,,  i  =  j  =  l,...7,  (2.65) 

Z 

SI  SJ 


Therefore,  denoting  the  Cartesian  velocity  components  of  the  “leader”  UAV  hy  X ,  Y  , 
and  Z ,  Eq.  (2.64)  is  rewritten  as 
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8=  -2 


A,e 

cl 


+  XlT^ss^ 


si 


+  . . .  + 


7  X 


\s7 


+ 


Yir^ee^ 


si 


s7 


(2.66) 


+  . . .  +  T^eSy  l^j  +  Z  [T^ee^  [i  +  ■  •  ■  +  ^7^^. 

and  thus  the  desired  Cartesian  velocities  that  the  UAVs  must  have  are  calculated  as  follows 

=  ~^x  +  •••  + 

¥“  ^-kjr^eej  +...  +  r,££,J  1  (2.67) 


Is7 


Y  \  1  Y 

Z‘=-kJr^Ee^ 


[^  +  ...  +  T^ee^ 


\s7 


s7 


The  error  gradient  at  the  leader  sensor  location  can  be  calculated  with  a  central  difference 
approach  using  the  measurements  of  the  followers 


(2.68) 


where  dl  is  the  distance  between  the  leader  UAV  and  a  follower,  assumed  to  be  equal  in  all 
the  directions  X ,  Y ,  and  Z .  The  error  gradients  at  locations  si  —  s3  and  s5  —  s7  are 
approximated  either  with  the  forward  difference  approach  or  with  the  backward  difference 
approach  using  the  neighbor  readings. 
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3  Numerical  Model 


In  this  chapter,  the  numerical  implementation  of  each  of  the  steps  in  the  approach  is 
discussed.  First,  the  chapter  provides  the  outline  of  the  whole  implementation  procedure. 
Second,  the  numerical  solution  of  the  advection-diffusion  equation  is  presented;  followed  by 
the  algorithm  for  the  eomputational  grid  adaptation.  At  the  end  of  the  ehapter,  the  details  for 
the  visualization  of  the  UAV  flight  are  given. 

3.1  Overall  Estimation-CFD  Implementation  Procedure 

The  mathematical  model  described  in  Chapter  2  has  been  implemented  numerieally  in  the 
3D  domain  as  shown  in  Figure  12,  with  winds  and  eddy  diffusivities  representative  of  the 
ambient  atmosphere  at  the  region  of  interest. 

The  UAV  starts  patrolling  the  domain  of  interest  following  a  eonstant-elimb  orbit  path 
with  predefined  desired  ground  speed  ,  course  angle  x'’’  flight  path  angle  Y  inputs 
to  the  dynamic  model  (2.42).  The  concentration  sensors  onboard  the  UAV  record  concentration 
(simulated)  data  according  to  its  spatial  location  0^  .  The  state  estimator  (2.18)  is  aetivated 

when  concentration  measurements  above  the  sensitivity  threshold  of  the  sensors  are  recorded. 


The  concentration  plume  is  estimated,  and  the  UAV  receives  new  estimator  performance-based 


Figure  12.  Real-time  gaseous  plume  estimation  approach  and  implementation. 
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control  signals  (2.58)  that  guide  the  vehicle  towards  the  areas  of  larger  estimation  error.  The 
transition  from  the  approach  to  the  mathematical  modeling  and  finally  to  the  implementation 
procedure  is  shown  in  Figure  12. 


3.2  Numerical  Solution  of  the  Advection-Diffusion  Equation 

The  process  model  advection-diffusion  equation  (2.13)  is  solved  numerically  in  order  to 
provide  numerical  data  in  the  absence  of  real  data  using  a  Finite  Volume  Method  (FVM)  [64], 
[66]  supplemented  with  the  Total  Variation  Diminishing  (TVD)  scheme.  The  estimator 
equation  (2.18)  is  solved  numerically  following  the  identical  FVM -TVD  procedure.  We  outline 
below  the  steps  in  3D  following  previous  2D  implementation  of  [24],  [25]. 

The  domain  under  consideration  Q  is  discretized  with  N  —  N ^  x  x  N ^  rectangular 

finite  volumes.  Integrating  Eq.  (2.13)  over  the  finite  volume  O, ,  with  surface  area  A  =  An 

tJK 

shown  in  Figure  13  yields 

f JJJ{c)d(l  +  §  F.dA  =  III  Sda  (3.1) 

%  ^  % 

The  vector  of  fluxes  F  represents  the  total  flux  through  the  surface  of  the  volume  and  is 
expressed  as  a  sum  of  advective  and  diffusive  fluxes: 


Figure  13.  A  finite  volume  at  grid  point  (i  j,k). 
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F  =  /"  +  /"  X  +  /"  +  /"  Y  +  //  +  /"  Z 


Based  on  Eq.  (2.13),  the  advective  and  diffusive  fluxes  in  X ,  Y  ,  and  Z  directions  are 
calculated  as  follows 

k'={cP.  fy=ic)V,  fy'={c)W 

fD  _  TZ  fD  _  TZ  fD  _  TZ 


•  f  =  -K 

1 


f^=—K  '  f  =  —K  '  f  =  —K  ' 

Jx  ^XX  ’  h  ^YY  QY  ’’  ^  QZ 

Eq.  (3.1)  is  replaced  by  its  discrete  form,  where  the  volume  integrals  are  expressed  as  the 

averaged  values  over  the  cell  and  where  the  surface  integral  is  replaced  by  a  sum  over  all  the 

bounding  faces  A^,  (  =  1, ...,  6  of  the  considered  volume 


■I  +  ^^Jk 


Due  to  the  rectilinear  (Cartesian)  grid,  the  calculation  of  the  dot  product  in  (3.4)  simplifies 


'  '  ijk  _ 

dt  Q 


F.lAt  -  FAA^^  +  F"A^^  -F^-A  .,  +  F  A  .  -  FA  ,  +  ^ 

ijk  ijk  ijk  ijk  ijk  ijk  ijk  ijk  ijk  ijk  ijk  ijk  I  i 


Using  Eq.  (3.2),  Eq.  (3.5)  can  be  rewritten  as 

ijk  ^ 

[fy  +  ^Y  )[^,  AJ  -  (//  +  fY  )|,^,  + 


^  )  ^,k  ^  ^ 


The  diffusive  flux  is  approximated  with  a  central  differencing  at  the  cell  interface  [66] 

nl^  \F  —  (c).  .  ,  ^\W  \W  (c).  .  ,  —  (c).  , 


=  -K 

ijk  ^ 


=  -K 

Jv  ..,  y 


=  -K 

ijk  2 


X.  , 

-X. 

hJ,k 

ijFl-,k 

Y 

i,j+l,k 

^hhk 

i,j,kFl 

i,j,k+l 

^hhk 

hj,k  fU  _  _T^ 

^  ijk  ^ 


=  -K 

’  Jy  Y 


=  -K 

5  Jz  ..,  ^^7 


^hhk 

^i—l,j,k 

-  F. .  , , 

hJ,k 

i,j—l,k 

i,j,k  z.  -  z. 

i,j,k  i,j,k—l 


The  advective  term  requires  special  attention  in  order  to  provide  monotonicity,  a  condition 
for  a  numerical  scheme  so  that  no  new  extrema  be  created  in  case  of  large  wind  speeds  near 
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large  gradients.  In  previous  work  [23]  the  upwind  scheme  was  implemented  in  2D.  In  the  3D 
implementation  of  this  dissertation  the  advective  flux  is  approximated  using  the  total  variation 
diminishing  (TVD)  scheme  [66],  [67].  The  algorithm  is  provided  below  for  a  simple  case 
where  V  =  W  =  {)  and  U  not  zero.  However,  the  created  code  assumes  the  possibility  of 
nonzero  wind  speeds  U,  V,  and  W.  The  wind  speeds  are  implemented  either  as  constants  or  as 
functions  of  special  variables  according  to  Eq.  (2.15). 


J  X 


Jx 


ijk 


W 


ijk 


>  0 


hj,k 


<0 


(3.8) 


>  0 


'U-  +-dr- 


'  ij,k  i,j,k  ■  2 


W 


hj,k 


<0 


The^(r)  is  the  limiter  function  of  the  local  ratio  of  upstream  to  downstream  gradient  r ;  the 
“+”  and  superscripts  refer  to  the  positive  or  negative  wind  direction  respectively  as  : 


r+  = 


r+  = 


^^h,j,k 

(^)7-l,i,A: 

.  ^i,j,k 

~  ^i-l,j,k  ^ 

,A:  ^^^7— 2,j,A: 
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{cl 
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\^hj,k 
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(^h,k  -  <c) 


7— l,j,A; 
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—  (c) 

\  /i+lj,k 

\  U.j.k 

X  ,  , 

-X,  ,, 

7+1, J, A: 

hJ,k 

X  ,  , 

-X  4 
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(c) 

—  (c) 

\  /i+lj,k 

\^u,j,k ; 

^  k  - 

X.  ^ 1 

7— l,j,A: 

(c). 

(c).  . , 

'  ' 7— l,j,A:  ^ 

(3.9) 


The  limiter  function  ^{rj  is  used  in  the  TVD  scheme  in  order  to  prevent  discontinuities  in  the 
obtained  solution.  When  ^  (r  j  =  1 ,  the  scheme  for  advection  term  reduces  to  the  average  of 

the  cell  centered  fluxes;  when  ^  (r  j  =  0 ,  the  scheme  reduces  to  the  first  order  upwind  “donor 

cell”  scheme.  For  the  proposed  estimation  algorithm,  the  limiter  function  is  chosen  as  the  Min- 
Mod  [68] 

lin  (r,  ij  if  r  >  0 


«(’■)  = 


mm( 
0 


if  r  <  0 

Assuming,  17  >  0 ,  and  substituting  Eq.  (3.7)  and  Eq.  (3.8)  into  Eq.  (3.6)  yields 


(3.10) 


29 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


-^  =  ^—\\{c).  .JJ.  ^-(c). -.U.  ,J- 
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The  limiter  function  ^  j  depends  on  the  value  of  concentration  (c)  and,  therefore,  varies 

with  time.  Consequently,  the  terms,  containing  multiplication  by  ^  j  can  be  considered  as  an 
additional  source  term,  which  is  referred  to  as  the  TVD  deferred  correction  source  term  and 
denoted  by  S’^'^ .  Thus,  Eq.  (3.11)  can  be  rearranged  to  yield 


-1 


-  -  -  \^/  •  I  1  -7  I  ^o\^/  •  1  -7  I  ^o\^/  •  ■  I  1  7  I  ^/l  \^/  •  •  17 
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(3.12) 


where  the  coefficients  a^,  i  =  1, . . . ,  6  are  calculated  based  on  Eq.  (3.1 1)  as  follows: 

E  (  IT  i 

K  K 

a,  = - "XmA - Ak^,  a  =  -U.  ^  - ^XhlA - 

1  V"  V  2  z—l.i.k  V"  V 

z+l,j,A:  i,j,k  i,j,k  t—lj,k 

V 

N  S 

3  Y  —  Y  ^  Y  —  Y  ^ 

i,jEl,k  i,j,k  i,j,k  i,j—l,k 

T  B 

^  _ _ hj,k  aT  n  — _ hhk  aB 

5  Y  _ Yi  6  Y  _ Yi 

i,j,k+l  i,j,k  i,j,k  ij,k—l 


(3.13) 


and  the  coefficient  a„  is  defined  as 


%  =  -«1  -  ^2  -  «3  -  -  «5  -  «6  +  (^,,jA,j,k  -  ^^-YjA,j,k 


(3.14) 
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The  deferred  correction  term  results  from  Eq.  (3.11) 


ijk 


s. 


DC 

ijk 


1  1 


I  2 


9  ^  ){^^'^i+l,j,k^i+l,j,k  ij.k)  ^i,j,k 


Substituting  Eq.  (3.7)  and  Eq.  (3.8)  into  Eq.  (3.6)  assuming  that  U  <  0  yields 
1 


dt 


n 


ij^k 


-K 


XX 


hj,k 


hj,k  X  —  X 

i+lj,k  iC,k 


A^.k- 

hJX 


^^\,j,k^i,j,k  ){^^\-l,j,k^i-l,j,k  ^^\,j,k^i,j,k)  ' 


-K 


XX 


i—l,j,k 


X  -X 

i,j,k  i—lj,k 


^  .  - 

hJX 


K 


YY 


N 


(c). 


i,j-\-l,k 


(c) 


i,j,k  aN  _\_  TC  I*^  ^^hjik  i.i— 


i,j,k  Y  —  Y 

i,j-\-l,k  iC,k 


-k  +  Kv 

i,j,k  YY 


Y  -  Y 

ij,k  ij—l,k 


ij-^,k  aS  _ 

hj,k 


K 


zz 


T  (c),,j,^+i  -  {c)i 


hhk  aT 


D,k  z.  —  Z. 

2j,A;+1  t,j,k 


A',,  +  K.. 


i,j,k  z. .—  z. , 

t,j,k  t,j,k—l 


i,j,k-l 

hj,k 


(3.15) 


(3.16) 


+  ‘S,. , 


hj,k 


Rearranging  Eq.  (3.16)  to  the  form  of  Eq.  (3.12)  results  in  the  following  coefficients 
i  —  1, . . . ,  6  and  : 


K 


U 


XX 


hjjk 


i+lC,k  Y"  _  Y 

2+1, j, A;  i,j,k 


K 


XX 


w 


hj:k 


A  a  = - 

i,j,k^  2 


hJ,k 


hj:k 


i—lj,k 


K, 


YY 


%=- 


K. 


iN 


YY 


hj:k  aS 

hj.k 


Ay...  a=  — 

Y  _  Y  ^  Y  —  Y 

i,jYl,k  i,j,k  i,j,k  2,j— 1,A; 

\T 


zz 


^ - 


i,j,k 


i,j,kYl 


ry  ^i,j,k'^  ^6 

- 


Kxz 

B 

hj:k 

^hj,k 

^i,j,k—l 

X.% 

hJ^k 


(3.17) 


%  S  S  ^4  ^5  %  [y iYlJ,k^iJ,k  ^ iC^ky.j.k 


(3.18) 


The  deferred  correction  term  S  ,  is  defined  as 

ijk 
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1  1 


1  i-k  I i,j,k  ^^\-\-l,j,k^  i-\-l,j,k)  ^i,j,k 


2^{''"w']{^^h-l,j,k^i-l,j,k  ^^h,j,k^i,j,k')  '^i,j,k 


(3.19) 


Combining  Eq.  (3.13)  and  Eq.  (3.17),  the  coefficients  a  and  a  are  written  in  the  general  case 


«i= 


«2=  ,_^,0) 


i-\-lJ,k  ij,k 


hj,k  j^E 

_  ^  ij,k 


i,j,k  t—lj,k 


(3.20) 


Combining  Eq.  (3.14)  and  (3.18),  the  coefficient  a  is  written  as 


a„  =  —a,  —  —  a— 

0  1  2  3  4  5  6 


or  (3-21) 


.„0J  +  -  ^max([/._^_.^,0j  +  min([/.^^^,0jj  A 

Einally,  combining  Eq.  (3.15)  and  Eq.  (3.19)  yields  the  general  expression  for  the  deferred 
correction  term  *5^*^ 


A  =  -jT  A“  -  (i  -  “.)« I  - 

ijk  Lv 

!(“.« (’'4  -  ('  -  “.)« J 


(3.22) 


where. 


a  =  1  for  (7.  ,  ,  >  0  and  a  =  1  for  [/.,  ,>  0 

e  z+lj,A:  w  i—l,j,k 

a  =  0  for  (7.  ,  .  <  0  and  a  =  0  for  [/  ,  ,  <  0 


(3.23) 


The  advection-diffusion  equation  (2.13)  (or  the  process  model)  is  supplemented  with  Neumann 
boundary  conditions,  which  implies  that  the  flux  at  each  physical  boundary  is  equal  to  some 
constant  value.  Consider  the  East  physical  boundary  i  —  N ^  of  the  domain  under 
consideration.  The  Neumann  boundary  condition  is  written  in  the  discretized  form  as  follows 


X)N^+\,i,k  X)N^,j,k  _  ^ 

X  —  X 

N  +lj,k  N  j,k 


(3.24) 
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where  E  are  given  values  for  all  j  =  1, . . . , N  and  k  —  . 

Expressing  {c)^  from  Eq.  (3.24)  and  substituting  it  into  Eq.  (3.11)  results  in  an 
additional  component  in  the  source  term 


K 


XX 


E,  At 


Vt 


N^,j,k 


(3.25) 


Substituting  {c] 


N^+l,j,k 


from  Eq.  (3.24)  into  Eq.  (3.16)  results  in  an  additional  component 


in  the  source  term 


J^k 


E, 


j,k 


K 


XX 


Ny,j,k 


-  X.  ,,  ,,  -X  ,,  17 


Nj.+l,j,k  N^,j,k  I  Nj.-\-l,j,k 


K 

N^,j,k 

Nx,J,k 


(3.26) 


Therefore  combining  Eq.  (3.25)  and  Eq.  (3.26)  yields  the  expression  for  an  additional 
component  in  the  source  term  taking  into  account  all  possible  wind  directions 


=  \E, 

J,k 


jA 


K 


XX 


NxJ,k 


-  X 


—  X 

^Nx+lJ,k  Nx,j,k 


mmiU 


Nx-\-l,j,k 


,0 


K 

Nx,j,k 

N  j,k 


(3.27) 


The  volume  + 1,  j,  j  lies  beyond  the  scope  of  the  domain  under  consideration  and 
its  coordinate  X  can  be  expressed  in  terms  of  X  as  shown  in  Figure  14.  The  value 
of  wind  speed  ^  can  be  found  by  linear  interpolation  using  the  values  of  the  wind  speed 


at  and 


(N 


• 

O 

• 

O 

• 

o 

•-H- 

-fko 

-^•=1,2, 


i—Nx  i—NyEl 


X 


Figure  14.  Finite  volume  discretization  at  the  boundary. 
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j  _  ^N^,j,k  ^N^-l,j,k 

N^+l,j,k  ~  ^  ^ 

N^,j,k  N^-l,j,k 


V  _  V  -\-TJ 

Ny+l,j,k  Ny,j,k]^  ^  Ny,j,k 


(3.28) 


The  term  S..^  in  Eq.  (3.12)  is  represented  by  a  vector  of  N  =  x  x  entries, 

which  are  all  zeros  except  for  the  single  volume  that  corresponds  to  the  source  location  0^  (t^ . 
Therefore  based  on  Eq.  (2.2), 


^  ^|0,  0,(t)^(ijk) 

hn{t),  e(t)e(ijk 


(3.29) 


The  numerical  solution  of  the  estimator  equation  (2.18)  follows  the  same  procedure  as  the 
advection-diffusion  equation  with  addition  of  the  output  injection  term  (2.19) 


'^jk  ~ 


0.  s,  (*)  ^  (ijk) 

>»r({c>(*.e)-(«>(‘.e)). 

(3.30) 


The  resulting  set  of  N  semidiscrete  equations  (3.12)  are  written  in  state  space  form  as 
follows 

x(f)  =  +  +  (3.31) 

where  the  source  term  is  expressed  as  the  product  of  the  vector  Q  (x,  E,  zj ,  which  represents 
the  location  information  of  the  source,  and  the  source  release  rate  uit).  The  vector 


X  =  (c)  =  (c)  , . . . ,  (c)  is  the  vector  of  states  that  represents  the  concentration  for  each  finite 


volume  in  the  computational  domain  (2 .  The  mapping  for  the  volumes  is  expressed  as  follows 

n  =  n  (i,  j,  k)  =  i  +  (^j- 1)  N^+(^k- 1)  N^N^, 

j  =  l,...,NY,  k  =  l,...,N^,  (3.32) 

\  =  i^).jk 


This  ordering  creates  matrices  P  and  P  from  Eq.  (3.31)  in  the  form 
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'^n  • 

: 

P 

TVDW 

..  P 

TVDIN 

Pn.  ■ 

PnN 

5 

P 

TVDnl 

..  P 

TVDnN 

Px.  ■ 

Pnn 

P 

TVDNl 

..  P 

TVDNN 

(3.33) 


The  matrix  P  is  a  7  diagonals  matrix  that  represents  the  finite  dimensional  advection 
diffusion  operator  with  entries  of  the  coefficients  a  from  Eq.  (3.17),  Eq.  (3.18),  Eq.  (3.20), 
and  Eq.  (3.21).  The  matrix  P^^,^  is  due  to  the  multiplication  by  the  limiter  function. 

The  system  of  ODEs  from  Eq.  (3.12)  and  Eq.  (3.31)  is  integrated  using  a  4*  order  Runge- 
Kutta  method  [64]  as 


X  = 

(0) 


1 


X  =  X - /\fP 

(1)  (0)  ^  RHS{0) 

^(2)  =  ^(0)  - 

X  ^  X  —  —  A  f  P 

(3)  (0)  2 

X  =  X  —  /\fP 

(4)  (0)  RHS{3) 

—  X 

(4) 


(3.34) 


The  Von  Neumann  method  [65]  is  used  for  stability  analysis  of  the  numerical  scheme.  In 
order  to  do  so,  we  consider  two  limiting  cases  for  the  limiter  function  ^  :  C  (^)  =  0  and 

(r  j  =  1 .  At  ^  (r  j  =  0 ,  i.e.  in  case  of  the  first  order  upwind  scheme  for  the  advective  flux. 


assuming  S..^^  =  0 ,  and  constant  eddy  diffusivities  and  wind  ( K 


XX 


^,J,k 


w 


=  K 


XX  ' 


YY 


hhk 


=  K 


YY 


YY  ’ 


K 


zz 


hhk 


=  K 


zz 


hhk 


=  A, 


zz  ■ 


and  (7,  ,  ^ 


V.  . ,  —W.  ^  ,  =  0 ),  Eq.  (3.1 1)  is  rewritten  as 

Z,7,rC  Z  —  i,7,rC  '*■ 
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d{c) 


ijk 


dt  n 


hj,k 


-(4_ 


■  1  , 

z— J 


+KA^.^ 

XX  i,j,k 


^KyyA^., 

YY  i,j,k 


iic) 

\  /z+l,j,A: 

^^\,j,k 

^^\,j,k 

X 

-X 

X 

-  X 

^  z-|-l,j,A; 

i,i,k 

i,j,k 

i-l,j,k  J 

iic) 

\  /z,j+l,A; 

Y 

-  Y 

Y 

-Y 

z,  j-|-l,fc 

hj^k 

^J,k 

i,j-l,k  ) 

f(c) 

\  /z,j,A:+l 

^^\,j,k 

z 

-  z 

z 

-Z 

^  z,  j,fc-|-l 

t,j,k 

t,j,k 

i,j,k-l  ) 

(3.35) 


In  the  uniform  grid,  such  that  X  , ,  —  X  ,  =  X  —  X  ,  ,  =  AX , 

^  z+lj,/c  ij^k  'f'^J^k  z— lj,/c 

K  -K  =y,  -K  =Ay,  Z,  -Z,  ,  ,  -Z  =  AZ,  and 

n  =AXAFAZ,  =AFAZ,  X""  =AXAZ,  =AXAr,  Eq.  (3.35)  is 

Ilk  ^  1.1  .k  ^  1.1  .k  ^  1.1  .k  ^  ±  Y  / 


i,j,k 

written  as 


dt 


AX 


+  X 


(^\+i,],k  ~  ^(^\,],k  +  (^\-i,],k 


XX 


I  +  K3~\k  I  ^  w/*  +  1 

^  YY  ^-^ZZ 


AX'  (3.36) 

('^)i.i.A:+l  “  2(C)  +  (C) 


AZ' 


We  also  represent  the  time  derivative  with  the  first-order  differencing 


ijk 


dt 


At 


(3.37) 


(7  At  ^  K^^At  ^  K^^At  ^  K^^At 

By  defining  as  the  CEL  number  a  = -  ,  f3  =  — ,  /5  =  — ,  /5  =  — 

AX  ^  AX'  AX'  ^  AX' 


and  using  Eq.  (3.37),  Eq.  (3.36)  is  rewritten  as 


+  ^r 


(3.38) 


Substituting  into  the  numerical  scheme  (3.38)  a  Eourier  mode  in  the  form 
(c)' where  ^  =  k^AX ,  i^^k^AY,  6  ^  k^AZ  ,  with  k^,  k^,  and  k^ 

being  the  wave  numbers  defined  as  k^  =  M^tt! ,  k^  =  M^tt! ,  andA:^  =  M^iTj , 
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where  ,  and  are  numbered  grid  points  in  the  X ,  Y ,  and  Z  directions 

accordingly,  yields 


ttI-I-I  Ii(p  Ij'ip  IkO  T  rl  liip  Ij'ifj  Ikd  (  t  rl  liip  Ij'ifj  Ikd 

V  e  e  e  =  V  e  e  e  —  (j\V  e  e  e  —  V  e  ^  ^  ^  ^ 


e  ‘'"e 


+13^  I  +  V‘e 


lj[i-l)^^lji;^lke 


(3.39) 


Simplifying  all  the  terms  in  Eq.  (3.39)  by  the  factor  and  determining  the 


yl+\ 


amplification  factor  as  (S'  =  yields 


G  =  l-a(l-  e-'^l  +  (3^  +  e'^A  +  X  +  e-^A  +  (3^  +  e'^A  - 

-  2{P^  +  0,- +  0,)  = 

1  —  (j(l  —  cos<^  +  / sin</7j  +  2(3^  cosip  +  2/?^  cosi/’  +  2(3^  cos9  — 

-  2(/?^  + 


For  If  =  'ip  =  9  =  7T , 


and  the  stability  condition 


becomes 


G  —  1  —  2cr  —  4  +  Py  Y  (3^ 


At  < 


(S'  <1 


1 


(3.41) 


(3.42) 


U 

AX 


+  2 


^XX  I  ^YY  I  ^ZZ 

AX'  AX'  AX' 


(3.43) 


For  cp  =  Ip  =  9  —  0 ,  G  =  1,  and  the  stability  condition  is  always  satisfied. 

At  ^  (rj  =  1 ,  i.e.  in  case  of  the  central  difference  scheme  for  the  advective  flux,  assuming 


=  0 ,  and  constant  eddy  diffusivities  and  wind  ( K 


XX 


hj,k 


K 


XX 


w 


hj:k 


K 


XX  ■ 
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K. 


YY 


\N 


YY 


S 


YY  ’ 


K 


zz 


T 


K 


zz 


B 


K, 


zz  ’ 


and 


^ ij,k  ^ ^  ’ 


V. . ,  =  W.  ^  . ,  =  0),  Eq.  (3.1 1)  is  written  as  follows 


i—lj,k 


d{c) 


ijk 


dt 


n. 


hj,k 


X  -  X,  , 

x+1  z— 1 


UA^-k 

hJ^k 


+k<xxAf,j,k 


+K^A^.. 

YY  t,j,k 


\ 


'  'i'd^k 


hj,k 


'  i—\,i,k 


z,j,A:  z— l,j,K 


Y  -Y 

i,j-\-l,k  i,j,k 


Y  -Y 

i,j,k  z,ji— l,fc 


(c).  —  (c).  (c).  —  (c).  . 

'  >i,j,k-\-l  '  'i.i.k  '  'i.T.k  '  'z.z 


z,j,A:+l 


'  hj,k 


i,j,k 


'  z,j,fc— 1 


hhk—\ 


(3.44) 


On  the  uniform  grid  it  is  written  as 


dt 


YK^ 


“  2(c)^  ^  +  (c)^_j 


'  i—lj,k 


+K^ 


2AX  '  ■ 

(*^)i.j  +  l.fc  ~  2(c),;.  ,-fc  +  (c),;.,_. 


hiik 


ij—l,k 


YY 


AF' 


+  K. 


(c). 


AX' 

iJ^kYl 


‘^(^)u,k  +  (^)Li,k-l 


ZZ 


AF' 


(3.45) 


Following  now  the  same  steps  as  before  in  Eq.  (3.37),  Eq.  (3.38),  Eq.  (3.39)  yields  the 
amplification  factor 

G  =  1  -  |(e^-  -  e-^")  +  +  e-^^)  +  p^  ^  - 

—  2  (^P^  +  p^  +  p^^  =  (3.46) 

1  —  al sin</7  +  2/3^  cos(f  +  2/3^  cosip  +  2/3^  cos0  —  2^/3^  +  Py  +  P^^ 


For  stability. 


G 


<1 


(3.47) 


Following  [65]  that  stability  condition  becomes 


<  2/3^ 

(^yY  (^,<- 


(3.48) 


which  provides  the  conditions  on  the  timestep. 
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At  < 


2K 

U‘ 


At  < 


‘  K  K  K  ^ 

XX  YY 


AX'  AX'  AZ' 


(3.49) 


The  application  of  the  TVD  numerical  scheme  is  verified  with  a  solution  for  the  “pure” 
advection  problem  Eq.(2.13)  with  V  =  W  =  0,  K —  0,  and  with  an 

instantaneous  line  source.  The  analytical  solution  for  this  case  is  given  by  the  step  function  and 
the  comparison  with  the  numerical  TVD  solution  with  the  upwind  scheme  is  shown  in  Figure 
15. 


3.3  Computational  Grid  Adaptation 

The  simulated  (concentration)  sensor  data  are  generated  by  solving  the  advection- 
diffusion  equation  (2.13)  on  a  very  fine,  uniform  grid  using  as  inputs  the  source  parameters 
(release  rate  and  position)  and  atmospheric  parameters.  Such  a  fine  uniform  grid  would  not  be 
applicable  for  the  desired  real-time  implementation  of  the  estimator  (2.18).  In  previous  work 
[19],  [24],  [25]  supported  via  FA9550-09-0469,  the  grid  was  adapted  by  a  priori  generating 
several  adapted  grids.  Each  predefined  grid  contained  an  area  of  a  refined  uniform  grid  that 
covered  25%  of  the  area  of  interest,  whereas  the  rest  of  the  domain  was  covered  with  a  coarse 
uniform  grid.  In  this  dissertation  in  order  to  reduce  the  computational  cost  for  the  estimator 

- ^ -  TVD 

Upwind 


Figure  15.  Verification  of  the  numerical  scheme. 
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simulation,  thus  obtaining  real-time  data,  the  computational  grid  is  adapted  dynamically 
throughout  the  process  [69].  The  adaptation  algorithm  uses  the  current  UAV  location  as  a 
center  of  a  refined  grid  with  the  grid  size  dl  in  all  directions.  Away  from  the  sensor,  the  grid 
nodes  are  following  a  parabolic  distribution.  The  resulting  three-dimensional  grid  is  shown  in 
Figure  16.  As  soon  as  the  UAV  crosses  the  distance  of  the  minimum  grid  size  dl  a  new  grid 
is  generated  and  the  numerical  solution  obtained  at  the  previous  time  step  is  interpolated  to  the 
new  grid  accordingly  using  trilinear  interpolation  method  [70].  The  total  number  of  nodes  for 
the  adapted  grid  remains  the  same  during  this  process. 

For  verification  of  the  proposed  grid  adaptation  algorithm,  the  numerical  solution  is 
compared  with  the  steady-state  analytical  solution  of  the  advection-diffusion  equation  for  a 

stationary  point  source  located  at  with  a  constant  release  rate  q  kg/s; 


dc 

The  steady  state  solution  i.e.  —  =  0  results  in  the  well-known  Gaussian  plume  formula 

dt 


which  for  a  stationary  source  is  given  as  [48],  [49] 


{c)(X,Y,1)  =  — - ^ ^ 


u 

V  CJ 

to 

h 

1 

to 

h 

1 

to 

2K 

XX 

V 

Kyy 

/ 

1 

to 

h 

1 

to 

h 

1 

to 

V 

^xx 

Kyy 

(X-X) 


K 


1/2 


(3.50) 


Figure  16.  Adapted  nonuniform  computational  grid. 
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The  steady-state  numerical  solution  is  achieved  at  time  t  ,  provided  that  the  system  is 

advection  dominated.  Numerical  and  analytical  solutions  are  compared  in  Figure  17:  to  the  left 
is  the  concentration  profile  obtained  on  the  uniform  grid  and  to  the  right  is  the  concentration 
profile  obtained  on  the  nonuniform  adapted  grid.  The  domain  under  consideration  Q,  has 
dimensions  4  km  x  4  km  x  1  km .  The  number  of  control  volumes  in  the  uniform  grid  is 
N  =  160  X 160  X  40 .  The  number  of  control  volumes  in  the  adapted  grid  is  =  40  x  40  x  10 
with  the  smallest  grid  size  h  .  =  12.5  m  and  the  effective  grid  size  (average  between  the 


maximum  and  minimum  grid  sizes)  =  100  m  .  The  results  are  shown  at  time  t  =  10  min. 

The  accuracy  of  the  numerical  method  is  based  on  the  spatial  error  norm  [71],  [72]  given 
by 


E  = 


EN 


numerical 


/\  Gaussian 
'  /  n 


n—1 


.c) 


Gaussian 


(3.51) 


The  log-log  plot  in  Figure  18  shows  the  relative  error  (3.51)  as  a  function  of  characteristic 
mesh  size.  Five  grid  levels  are  considered  for  the  uniform  grid  (shown  in  blue  in  Figure  18), 
with  h  =  200,  100,  50,  25,  and  12.5  m.  The  effective  grid  size  for  nonuniform  grid  (shown 

in  red  in  Figure  18)  ranges  from  25  m  up  to  200  m,  with  the  smallest  grid  size  h  .  =  12.5  m. 


The  analysis  shows  higher  accuracy  of  the  non-uniform  adapted  grid  at  larger  characteristic 

- -  Numerical  solution 

•  Analytical  solution  (Gaussian  plume) 


Figure  17.  Grid  adaptation  verification:  numerical  solution  vs.  analytical  formula. 
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Figure  18.  Relative  error  as  a  function  of  characteristic  mesh  size. 


mesh  sizes.  The  larger  is  the  characteristic  mesh  size,  the  smaller  is  the  number  of  control 
volumes  N  to  be  used.  Smaller  N  results  in  the  faster  implementation  of  the  numerical 
scheme.  Therefore,  the  use  of  the  grid  adaptation  is  beneficial  and  maintains  the  required 
accuracy  of  the  estimation  provided  that  the  number  of  control  volumes  is  limited  due  to  the 
necessity  to  implement  on-line  computations. 


3.4  Numerical  Simulation  of  UAV  Dynamics  and  Guidance 

The  UAV  state  vector  follows  Eq.  (2.42)  and  is 


T 


(3.52) 


The  commanded  input  vector  follows  Eq.  (2.58)  and  is 


T 


(j)  T 


(3.53) 


In  order  to  integrate  in  time  numerically  the  state  vector,  given  the  vector  of  commands, 
it  is  necessary  to  calculate  all  the  values  used  in  the  dynamic  model  (2.42).  Given  the  values 
for  the  ground  speed  U  ,  course  angle  %  and  flight  path  angle  7  at  the  current  time  step,  the 

airspeed  magnitude  U  can  be  expressed  from  Eq.  (2.33).  The  heading  angle  'ip  is  then 

calculated  using  Eq.  (2.34)  and  Eq.  (2.35).  Finally,  the  system  of  ODEs  (2.42),  that  can  be 
written  in  the  vector  form  as 
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i>  =  f{t),  (3.54) 

supplemented  with  Eq.  (2.43)  and  Eq.  (2.44)  is  integrated  using  4*  order  Runge-Kutta  method 
to  yield  the  new  state  veetor: 

Pk+i  ^  Pk  +  ~^(K  +  +  2^3  +  ^4) 

where 


K  =  f{h) 


K  =  f 


K=f 


.  At 

t  H - 

k  2 

.  At 
t,  H - 

k  2 


^4  =/k  +  A^) 


(3.56) 


For  the  flight  visualization,  the  Euler  angles  9  and  ip  should  be  calculated  in  addition  to 
bank  angle  (p ,  which  is  a  given  control  input.  To  calculate  the  heading  angle  ip ,  the  same  steps 
are  implemented  as  before,  i.e.  the  airspeed  magnitude  F  and  the  air-mass-referenced  flight 
path  angle  7^  are  recalculated  using  Eq.  (2.33)  and  Eq.  (2.34)  accordingly  to  be  used  in  Eq. 

(2.35)  The  pitch  angle  6  results  from  Eq.  (2.30),  where  the  angle  of  attack  a  is  calculated 
based  on  the  control  input  C and  Eq.  (2.45). 
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4  Summary  of  Simulation  Results 


4.1  Simulation  Parameters 


The  approach  presented  is  applied  to  gas  releases  in  the  atmosphere  under  conditions 
specified  in  Table  1  for  a  stationary  source  and  a  moving  aerial  source. 

The  UAV  is  modeled  using  physical  parameters  for  the  Aerosonde®  UAV  shown  in  Table 
2  [45]  and  summarized  below. 

The  domain  under  consideration  has  dimensions  Q  =  4km  x  4km  x  1km .  The 
advection-diffusion  equation  (2.13)  is  discretized  with  N  =  200  x  200  x  50  finite  volumes. 
The  number  of  finite  volumes  for  the  discretization  of  the  estimator  (2.18)  is 
=  80  X  80  X  20  with  the  minimum  grid  size  dl  =  20  m.  The  sensor  is  assumed  to  provide 
noiseless  measurements  with  the  time  response  of  1.5-2  seconds  and  concentration  range 


10 


-10 


kg/m^ 


Table  1.  Atmospheric  parameters 


Parameter 

Value 

Parameter 

Value 

Density  p 

1.2922  kg/m^ 

Eddy  diffusivity  Kxx 

20  mVs 

Wind  speed  U 

[7,9]  m/s 

Eddy  diffusivity  Kyy 

20  mVs 

Wind  speed  V 

0  m/s 

Eddy  diffusivity  Kzz 

10  mVs 

Wind  speed  W 

0  m/s 

Table  2  UAV  Specifications. 

Parameter 

Value 

Parameter 

Value 

Mass  M 

13.5  kg 

Oswald  efficiency  factor  e 

0.9 

Planform  area  S 

0.55  m^ 

Parasitic  drag  coefficient  c  ^ 

0.0437 

Wingspan  b 

2.8956  m 

Lift  coefficient  in  stall 

1.632 

Cruise  speed 

25-30  m/s 

Maximum  thrust  Tmax 

50  N 

Maximum  speed 

33  m/s 

Maximum  bank 

30P 
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4.2  Aerosonde®  UAV 


The  Aerosonde  UAV  shown  in  Figure  19  has  a  wingspan  of  2.9  m  and  weight  of  13-15 
kg.  It  is  powered  by  24cc  fuel  injected,  premium  unleaded  gasoline.  The  cruise  speed  ranges 
from  20  to  40  m/s.  The  maximum  flight  duration  is  about  30  hours  with  the  distance 
approximately  3,000  km.  The  maximum  altitude  ranges  from  0.1  to  6  km  depending  on 
payload.  The  payload  is  up  to  2  kg.  The  Aerosonde  UAV  is  navigated  by  GPS  and 
communicates  via  UHF  radio  or  LEO  satellite.  One  distinguishing  feature  of  the  Aerosonde 
UAV  is  the  rear  propeller,  which  allows  for  atmospheric  measurements  before  air  is  disturbed 
by  propeller. 

The  advantages  of  the  Aerosonde  UAV  for  our  application  are  the  following:  it  was 
designed  specifically  with  scientific  research  applications  in  mind:  data  is  high  in  quality  as  a 
result  and  all  parameters  related  to  the  UAV  are  under  control  by  the  user;  the  relatively  slow 
speed  of  the  Aerosonde  (20  to  40  m/s)  allows  instruments  onboard  to  collect  data  at  a  great 
sampling  rate;  it  is  able  to  fly  very  close  to  the  surface  thus  collecting  very  high  resolution 
data. 

Among  the  disadvantages  the  most  crucial  are  the  following:  the  total  payload  when  fully 
fueled  is  only  2  kg,  limiting  the  weight  of  instruments  the  Aerosonde  can  carry;  the  Aerosonde 
does  not  have  the  capability  to  detect  other  UAVs,  aircrafts,  or  other  obstructions  in  order  to 
avoid  them. 


Figure  19.  Aerosonde®  UAV  [45].  Copyright  2002,  John  Maurer  iniaurer@hawaii.edu 
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The  Aerosonde®  UAV  finds  many  applications  such  as  validation  of  satellite-derived 
products:  e.g.  passive  microwave  sea  ice  concentrations  and  sea  surface  temperatures;  support 
of  in  situ  data  collection  efforts;  imaging  of  sea  ice  in  polar  regions  where  conditions  are 
hazardous;  studying  conditions  leading  to  icing  of  an  aircraft  body  during  ~0°C  temperatures 
and  high  humidity;  development  of  technologies  to  prevent  these  conditions  on  aircraft 
carrying  and  being  flown  by  humans  (icing  can  potentially  cause  an  airplane  to  stop 
functioning  while  in  flight);  atmospheric  profiling  and  flux  studies;  mapping  sea  surface 
temperature  at  high  resolutions  using  an  infrared  thermometer.  In  addition,  micro-SAR  can  see 
through  clouds  and  in  the  dark,  thus  being  a  good  tool  for  mapping  sea  ice  in  polar  regions  at 
very  a  high  resolution  (~l-2  meters).  The  laser  altimeter  is  used  for  producing  highly  detailed 
digital  elevation  models,  which  can  be  applied,  for  example,  to  study  the  mass  balance  of 
Greenland.  Meteorological  data  can  be  collected  for  weather  prediction  and  weather  studies. 
Heavy  storms  can  be  flown  into  for  these  purposes  as  well,  even  into  hurricanes. 

4.3  Single  Sensor 
4.3.1  Stationary  Source 

A  stationary  gas  source  is  placed  inside  the  domain  at  a  point  with  coordinates 
0^  (fj  =  0^  =  (l.2km, 2km,0.5kmj  releasing  material  continuously  with  u  —  1  kg/s.  The 


- Source  location 

- UAV  trajectory 


Figure  20.  UAV  trajectory  for  stationary  source. 
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UAV,  having  no  information  about  the  source  location  and  strength,  starts  patrolling  the 
domain  at  a  point  with  coordinates  (oj  =  (2.8km,  0.8  km,  0.2  kmj ,  following  a  constant- 

climb  orbit  path  with  flight  path  angle  7  =  8°  and  orbit  radius  equal  to  0.4  km.  After 
approximately  280  seconds  of  simulation,  the  sensor  detects  a  nonzero  concentration,  which 
results  in  the  activation  of  the  estimation  algorithm  and  thus  switches  the  control  of  the  UAV. 
Since  a  nonzero  concentration  is  detected,  the  UAV  receives  the  control  signals  resulted  from 
the  desired  ground  speed  components  (Eq.  (2.58)),  which  guide  the  UAV  towards  areas  of 
larger  estimation  error.  In  the  case  of  continuous  gas  release  with  constant  release  rate,  this 
area  contains  the  source  location.  Figure  20  illustrates  the  resulting  UAV  path  in  horizontal 
plane,  vertical  plane,  and  in  three-dimensional  view.  Because  of  the  physical  constraints 
imposed  on  the  UAV  motion  and  the  presence  of  a  strong  wind  from  west  to  east,  the  UAV 
may  leave  the  source's  proximity  and  detect  zero  concentration.  If  the  sensor  reads  zero 
measurements  after  the  estimator  activation,  the  UAV  is  driven  towards  the  point  of  maximum 
estimated  concentration  provided  by  Eq.  (2.18).  Figure  5  shows  the  resulting  trajectory 
adjustment  due  to  this  condition  in  the  source's  proximity. 

The  performance  of  the  estimator  is  examined  by  calculating  the  RMS  error  of  the 
estimated  concentration  vector.  Although  the  source  localization  is  not  the  primary  goal  of  the 
present  work,  the  fact  that  the  guidance  scheme  moves  the  UAV  towards  the  location  of  the 


Figure  21.  Estimator  performance  analysis  for  a  stationary  source.  (Left)  the  RMS  error  vs.  time; 
(right)  distance  between  UAV  and  source. 
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Figure  22.  Comparison  between  estimated  plume  concentration  (left)  and  simulated  concentration 
(right)  for  a  stationary  aerial  source  at  various  times. 

continuously-releasing  source  proves  that  the  control  signals  indeed  guide  the  sensor  towards 
the  area  of  larger  estimation  error.  Therefore,  every  time  the  UAV  sensor  takes  a  measurement, 
the  estimator  provides  a  more  accurate  reconstruction  of  the  gaseous  plume.  The  RMS 
concentration  error  is  shown  in  Figure  21  (left)  as  a  function  of  time  and  the  evolution  of 
distanee  between  the  souree  and  the  UAV  is  shown  in  Figure  21  (right).  Oscillations  in  the 
RMS  error  are  due  to  the  time  that  the  UAV  takes  to  reach  the  next  desired  location  and  the 
sensor  time  response.  Spikes  at  the  end  of  both  plots  in  Figure  21  are  due  to  the  UAV  looping 
in  the  source's  proximity. 
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7E-07 


Figure  24.  Estimator  performance  for  a  moving;  (Left)  RMS  error  vs.  time;  (Right)  distance 

between  the  UAV  and  the  source. 


The  estimated  concentration  of  the  gaseous  plume  is  compared  with  the  simulated  data 
and  presented  in  Figure  22  for  different  time  instances.  It  is  seen  that  by  the  end  of  the 
simulation  time  the  shape  of  the  estimated  plume  becomes  closer  to  the  shape  of  the  “real” 
plume. 


Figure  23.  UAV  trajectory  for  a  source  moving  along  an  arc  trajectory. 
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4.3.2  Source  Moving  along  an  Arc  Trajectory 

In  this  application,  the  gas  source  is  moving  along  an  arc  trajectory  maintaining  a  constant 
altitude  with  the  rate  of  release  u  =  1  kg/s.  The  UAV  starts  patrolling  the  domain  at  a  point 

with  coordinates  0^  =  ^3.2km,  2.0  km,  0.2  km j  following  the  constant-climb  orbit  path 

with  the  flight  path  angle  7  =  3°  and  the  orbit  radius  equal  to  0.4  km.  The  activation  of  the 


Figure  25.  Comparison  between  estimated  plume  concentration  (left)  and  simulated  concentration 
(right)  for  a  source  moving  along  an  arc  trajectory. 
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estimator  occurs  after  190  seconds  of  simulation,  when  the  sensor  detects  a  nonzero 
concentration.  The  total  simulation  time  is  420  seconds. 

Figure  23  shows  the  resulting  UAV  trajectory.  It  is  seen  that  the  UAV  is  able  to  reproduce 
approximately  the  source's  trajectory.  However,  it  stays  at  some  distance  to  the  east  of  the 
source,  which  is  explained  by  the  wind  direction.  The  evolution  of  the  RMS  concentration  error 
and  the  distance  between  the  UAV  and  the  source  are  shown  in  Figure  24.  In  this  case,  the 
UAV  follows  the  source  and  after  detecting  a  concentration  with  its  sensor,  the  distance 
between  the  UAV  and  the  source  decreases  until  the  end  of  the  simulation  process. 

Figure  25  compares  the  estimated  concentration  of  the  gaseous  plume  to  the  simulated 
data  (plant)  at  different  time  instances.  Once  again,  it  is  observed  that  by  the  end  of  the 
simulation  time  the  shape  of  the  estimated  plume  becomes  closer  to  the  shape  of  the  “real” 
plume. 


4.3.3  Source  Moving  across  the  Domain  with  Altering  Altitude 

In  this  application,  the  gas  source  is  moving  across  the  domain,  with  its  altitude  changing 
as  a  sine  function.  The  rate  of  release  u  =  1  kg/s.  Unlike  previously  considered  examples 
(where  the  wind  speed  was  taken  to  be  7  m/s),  the  wind  speed  for  this  case  is  9  m/s.  The  UAV 
starts  patrolling  the  domain  at  a  point  with  coordinates  0^  (o)  =  ^2.4km,2.8km,0.2kmj 


Figure  26.  UAV  trajectory  for  a  source  moving  across  the  domain. 
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Figure  27.  Estimator  performance  analysis  for  a  source  moving  across  the  domain.  (Left)  the  RMS 
error  vs.  time;  (right)  distance  between  UAV  and  source. 

following  the  constant-climb  orbit  path  with  the  flight  path  angle  7  =  3°  and  the  orbit  radius 

equal  to  0.4  km.  It  takes  approximately  140  seconds  for  the  sensor  to  detect  a  concentration 

above  the  threshold.  The  total  simulation  time  in  this  case  is  480  seconds. 

Figure  26  shows  the  resulting  UAV  trajectory.  The  UAV  is  able  to  follow  the  source. 
However,  at  the  end  of  the  simulation  time  it  tends  to  move  upwind  (recall  that  the  wind  speed 
is  9  m/s,  whereas  the  Aerosonde®  UAV  cruise  speed  is  about  30  m/s).  The  estimator 
performance  analysis  is  presented  in  Figure  27.  As  time  increases,  the  RMS  error  decreases, 
which  is  indicative  of  convergence  of  the  proposed  approach.  For  a  couple  of  moments  during 
the  simulation  time,  the  distance  between  the  sensor  and  the  source  comes  to  zero.  This  requires 
additional  algorithms  for  collision  avoidance. 

Figure  28  presents  visualization  of  the  “real”  and  estimated  plumes  for  several  time 
moments. 
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Figure  28.  Comparison  between  estimated  plume  concentration  (left)  and  simulated  concentration 
(right)  for  a  source  moving  across  the  domain. 

4.3.4  Ground-Based  Source 

In  this  application  the  real-time  prediction  of  gas  contaminant  concentration  from  a  ground 
intruder  using  UAV  is  examined.  In  this  ease  the  eomputational  domain  considered  has 
dimensions  different  from  the  other  tested  seenarios  and  equal  to  3000m  x  3000m  x  600m . 
The  wind  speed  is  ealculated  aceording  to  Eq.  (2.15)  with  Z^  —  ^m  and  m  =  0.15; 

V  =  W  =  0.  The  estimator  is  approximated  with  =  30  x  30  x  40  volumes,  so  that  the 
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computational  grid  is  adapted  only  in  the  X  and  Y  directions  and  remains  uniform  in  the  Z 
direction.  The  minimum  grid  size  in  the  adapted  grid  is  c?/  =  15  m .  The  sensor  is  assumed  to 
provide  noiseless  measurements  with  1  second  time  response. 

In  the  examined  scenario,  the  intruder  is  moving  across  the  domain  from  the  North-West 
to  the  South-East  comer,  releasing  gas  with  a  rate  of  1  kg/s.  The  UAV  starts  patrolling  the 
domain  downwind  following  a  circular  path  of  400  m  radius,  at  an  altitude  of  120  m.  It  takes 
approximately  190  seconds  for  the  plume  to  reach  the  sensor.  After  that,  the  UAV  receives  the 
control  signals  to  follow  the  desired  Cartesian  velocities  (2.55).  The  UAV  can  fly  as  low  as  10 
m.  If  at  any  moment  after  the  activation  of  the  estimator,  the  sensor  has  a  zero  reading,  it  means 
that  the  UAV  moved  upwind  and  lost  the  plume.  In  this  case,  the  guidance  scheme  is  adjusted 
so  that  the  sensor  follows  the  maximum  estimated  concentration,  thus  returning  to  an  area  of 
nonzero  concentration  measurements.  Figure  31  shows  the  resulting  UAV  trajectory.  As  it  is 
seen,  the  UAV  is  moving  towards  the  intmder  with  oscillations  caused  by  the  wind  direction, 
the  sensor  time  response,  and  the  time  required  for  the  UAV  to  complete  a  maneuver.  The 
loops  in  the  trajectory  are  due  to  the  fact  that  the  sensor  leaves  the  plume  moving  upwind  and 
it  takes  several  time  iterations  to  return  it  back  to  the  source’s  proximity. 
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Figure  29.  UAV  trajectory:  red  line  -  UAV  path;  green  line  -  intruder. 
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Figure  30.  Distance  between  the  sensor  and  the  source  vs.  time. 

The  estimator  performance  is  quantified  by  the  behavior  of  the  normalized  RMS  error, 
which  measures  the  difference  between  the  estimated  and  true  concentration  field,  and  the 
distance  between  the  sensor  and  the  intruder.  Figure  3 1  shows  the  normalized  RMS  error  as  a 
function  of  time.  As  it  is  observed,  the  RMS  drops  rapidly  after  the  estimator  activation  and 


Figure  31.  Normalized  root-mean-square  error  vs.  time. 

55 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Concentra(ion (kg/m'):  IE- 10  6E-I0  2E-09  8E-09  4E-08  IE-07 


Figure  32.  Evolution  of  the  estimated  plume  from  240s  to  480s  with  80s  increments. 

stays  on  average  within  10%.  The  oscillations  are  due  to  the  fact  that  the  source  is  releasing 
material  continuously.  The  spikes  around  360s,  400s  and  460s  are  explained  by  the  fact  that 
the  UAV  loses  the  source’s  proximity,  and  it  takes  several  seconds  for  it  to  relocate  the  plume. 
Figure  30  shows  the  distance  between  the  sensor  and  the  source  as  a  function  of  time.  It  is  seen 
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that  the  sensor  approaches  the  gaseous  source,  which  proves  that  the  proposed  guidance 
scheme  indeed  moves  the  UAV  towards  the  larger  estimation  error  that  coincides  with  the 
source  location  in  the  case  of  continuous  release. 

Figure  32  shows  the  real-time  estimates  of  gas  concentration  from  a  ground  mobile 
source. 


4.4  Multiple  Sensors 

For  the  “practical”  application  of  the  suggested  estimation  approach  the  concentration 
gradients  are  estimated  using  multiple  sensor  readings.  In  Section  2.4.6  we  suggest  to  take 
advantage  of  the  additional  measurements  and  to  include  them  into  the  estimator  (2.18).  The 
performance  of  the  estimation  algorithm  is  tested  numerically  for  this  case.  The  sensors  are 
assumed  to  be  attached  to  different  UAVs  that  maintain  a  rigid  formation  during  the  whole 
simulation  time.  The  distance  dl  between  the  leader  UAV  and  a  follower  is  equal  in  all  the 
directions  X ,  Y  ,  and  Z ,  and  is  dictated  by  the  numerical  resolution  of  the  process  model 
equation  (2.13).  However,  in  a  real  life  situation,  such  parameter  would  be  chosen  based  on 

the  advection  length-scale  —Ut  and  diffusion  length-scale  —  4^  in  order  to  give  the 
meaningful  gradient. 


Figure  33.  Comparison  of  the  trajectory  in  case  of  the  concentration  estimation  with  a  single  sensor 

and  with  multiple  sensors. 
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Figure  34.  The  normalized  concentration  error  versus  time  for  the  case  of  a  single  sensor  and  the 

case  of  multiple  sensors. 

In  the  first  considered  case  the  source  is  moving  across  the  domain.  The  wind  speed  is  9 
m/s,  the  total  simulation  time  is  360  seconds.  The  results  are  shown  in  comparison  with  the 
single  sensor  case,  provided  that  a  single  sensor  has  full  knowledge  of  concentration  and 
concentration  gradients  at  its  location. 

Figure  33  shows  the  resulting  trajectory  for  a  single  sensor  and  the  sensor  formation.  In 
the  first  case,  the  trajectory  results  from  the  guidance  scheme  (2.55),  and  in  the  second  case  it 
results  from  Eq.  (2.67)  with  constant  F  ’s.  It  is  seen  that  the  flying  formation  is  more  successful 
in  chasing  the  source  by  the  end  of  the  simulation  time.  The  corresponding  distance  between 
the  sensor  and  the  source  is  shown  in  Figure  35. 

The  concentration  error  normalized  by  the  total  amount  of  material  released  is  shown  in 
Figure  34.  Specifically,  the  last  minute  of  the  simulation  is  zoomed  in  order  to  clearly  show 
the  oscillations.  It  is  seen  that  the  oscillations  are  smaller  in  case  of  multiple  sensors. 

Another  example  is  considered  for  application  of  the  multiple  sensors.  In  this  case,  the 
source  is  moving  along  an  arc  trajectory,  similar  to  Section  4.3.2.  The  simulation  time  in  this 
case  is  420  seconds.  The  trajectory  and  the  corresponding  distance  between  the  sensor  and 
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Single  sensor 
Multiple  sensors 


Figure  35.  The  distance  between  the  sensor  and  source  versus  time  for  the  case  of  a  single  sensor 

and  the  case  of  multiple  sensors. 

source  for  this  case  is  shown  in  Figure  37  and  Figure  36  accordingly.  Similar  to  the  previous 
example,  the  flight  formation  improves  the  estimation  approach  by  the  end  of  the  simulation 
time  as  compared  to  a  single  sensor. 
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Figure  36.  The  distance  between  the  sensor  and  source  versus  time  for  the  case  of  a  single  sensor 
and  the  case  of  multiple  sensors;  the  source  is  moving  along  an  arc  trajectory. 
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Figure  37.  Comparison  of  the  trajectory  in  case  of  the  concentration  estimation  with  a  single  sensor 
and  with  multiple  sensors;  the  source  is  moving  along  an  arc  trajectory. 
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5.  Overview  and  Conclusions 


5.1  Overview 

This  work  expanded  the  earlier  work  [18],  [19],  [20],  [21],  [22],  [23],  [24],  [25]  supported 
by  FA9550-09-0469  on  a  coupled  CFD-control  approach  that  provides  the  real-time  estimate 
of  a  gaseous  plume  in  the  atmosphere  presented  by  an  UAV  equipped  with  sensors.  These 
developments  involve  the  mathematical  model,  the  computational  model  and  the  applications, 
and  have  been  presented  also  in  [39],  [40],  [41],  [42],  The  mathematical  model  for  the  process 
of  gas  release  was  extended  to  a  3D  unsteady  advection-diffusion  equation  and  it  was 
implemented  numerically  using  a  finite  volume  method.  The  numerical  implementation  was 
supplemented  with  a  total  variation  diminishing  scheme  for  the  advection  term.  The  numerical 
code  created  in  Fortran  90  considers  different  source  trajectories  and  release  types  (pulsed  or 
continuous),  and  different  models  for  winds  and  eddy  diffusivities.  The  sensor  model  was 
expanded  to  a  3D  case  and  implemented  given  a  specific  response  time  and  sensitivity  range. 
A  3D  dynamical  model  for  the  UAV  was  developed  and  implemented  in  the  code.  The 
guidance  law  based  on  the  Lyapunov  redesign  methods  was  modified  in  such  a  way,  that  the 
Lyapunov  function  resulted  in  the  desired  direction  for  the  vehicle.  A  lower-level  controller  in 
turn  provided  the  UAV  with  the  inputs  that  a  pilot  commonly  controls:  thrust,  lift,  and  bank 
angle.  For  the  purpose  of  numerical  efficiency  in  the  real-time  tracking  problem,  the  algorithm 
for  the  computational  grid  adaptation  was  developed.  According  to  this  algorithm,  the 
computational  grid  is  adapted  continuously,  thus  keeping  an  area  of  high  resolution  near  the 
sensor.  The  desired  direction,  namely  Cartesian  velocity  components  for  the  UAV,  was 
expressed  in  terms  of  the  concentration  error  gradient.  The  concentration  error  gradient  was 
calculated  using  multiple  concentration  sensors,  which  were  assumed  to  be  attached  to 
different  vehicles  that  maintain  a  rigid  flying  formation  around  a  leader-UAV.  The 
measurements  from  additional  sensors  were  included  in  a  3D  estimator  model  via  the  use  of 
the  observation  operator  given  by  the  7-dimensional  vector.  The  approach  was  tested  as  applied 
to  different  source  trajectories  using  physical  parameters  of  the  Aerosonde®  UAV.  The 
applications  included  detection  of  a  ground-based  intruder,  and  an  aerial  source.  In  addition, 
the  case  of  a  flying  formation  for  estimation  of  the  concentration  plume  as  well  as  the  source 
proximity  was  examined. 
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Figure  38  shows  the  main  steps  of  the  approach  implementation.  The  plume  is  assumed  to 
be  generated  by  a  moving  aerial  source  with  unknown  strength  and  location,  and  is  modeled 
by  the  unsteady  advection-diffusion  equation  with  winds  and  eddy  diffusivities.  The  UAV 
equipped  with  concentration  sensors  patrols  the  domain  and  collects  sensor  readings.  As  soon 
as  the  sensor  reads  the  data  above  a  certain  threshold,  the  estimator  is  activated.  The  estimated 
concentration  results  from  the  numerical  solution  of  the  advection-diffusion  equation 
supplemented  with  the  concentration  sensor  readings.  The  approach  couples  the  UAV 
guidance  with  the  performance  of  the  estimator  through  the  state  (concentration)  estimation 
error.  Via  the  appropriate  choice  of  the  Lyapunov  function  for  the  resulting  state  estimation 
error  the  estimation  scheme  provides  the  desired  Cartesian  velocities  of  the  UAV.  A  lower- 
level  controller  in  turn  accepts  these  desired  velocities  as  reference  inputs  and  provides  the 
control  signals  to  the  UAV  in  order  to  follow  the  desired  velocities.  As  a  result,  the  UAV 
moves  towards  the  areas  with  the  maximum  estimation  error.  The  finite-volume  discretization 
of  the  estimator  incorporates  a  second-order  total  variation  diminishing  scheme  (TVD)  for  the 
advection  term.  To  achieve  greater  computational  efficiency,  the  computational  grid  for  the 
estimator  is  adapted  dynamically  with  local  grid-refinement  centered  at  the  UAV  location. 


Figure  38.Flow  chart  of  the  approach. 
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5.2  Conclusions 

The  mathematical  model  for  the  advection-diffusion  equation,  for  the  point-sensor 
measurements,  as  well  as  for  the  Luenberger  observer  estimator  was  expanded  to  a  3D  case. 
The  3D  dynamical  model  was  developed  for  the  UAV.  The  guidance  law  for  the  UAV  was 
modified  in  order  to  provide  the  control  inputs  in  terms  of  the  values  that  a  pilot  commonly 
controls:  engine  thrust,  lift,  and  bank  angle.  The  gradient  measurement  issue  was  addressed  by 
using  multiple  sensors.  The  estimator  was  modified  to  include  multiple  sensor  readings. 
Numerically,  the  discretization  of  the  advection-diffusion  equation  was  improved  by  adding 
total  variation  diminishing  scheme  for  the  advection  term.  In  addition,  the  equation  was  solved 
using  more  sophisticated  computational  grid  adaptation  approach.  The  computational  grid  was 
adapted  continuously  during  the  simulation  time  to  provide  higher  resolution  near  the  sensor 
location.  The  3D  UAV  dynamical  model  was  implemented  numerically,  and  the  aircraft 
response  to  the  control  inputs  based  on  the  estimator  performance  was  visualized.  The 
approach  was  tested  using  physical  parameters  of  the  Aerosonde®  UAV  for  different  source 
trajectories  and  atmospheric  conditions. 
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